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The strangeness contribution to the electric and magnetic properties of the nu- 
cleon has been under investigation experimentally for many years. Lattice Quantum 
Chromodynamics (LQCD) gives theoretical predictions of these measurements by 
implementing the continuum gauge theory on a discrete, mathematical Euclidean 
space-time lattice which provides a cutoff removing the ultra-violet divergences. In 
this dissertation we will discuss effective methods using LQCD that will lead to 
a better determination of the strangeness contribution to the nucleon properties. 
Strangeness calculations are demanding technically and computationally. Sophisti- 
cated techniques are required to carry them to completion. In this thesis, new theo- 
retical and computational methods for this calculation such as twisted mass fermions, 
perturbative subtraction, and General Minimal Residual (GMRES) techniques which 
have proven useful in the determination of these form factors will be investigated. Nu- 
merical results of the scalar form factor using these techniques are presented. These 
results give validation to these methods in future calculations of the strange quark 
contribution to the electric and magnetic form factors. 
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CHAPTER ONE 
Introduction 

In physics today there exist four fundamental forces: the strong force, the weak 
force, electromagnetism and gravity. The focus of this thesis is the strong force 
and related particles. Quantum Chromodynamics (QCD) is the study of the strong 
interaction. 

A hadron is a particle constructed of quarks and gluons, which are the fun- 
damental strong force particles. There exist six different flavors of quarks: up (u), 
down (d), strange (s), charmed (c), bottom (b), and top (t). Of these six quarks 
the up, down, and strange are known as the light quarks. The up and down quarks 
have masses of a few MeV while the strange quark has a mass of approximately 120 
MeV. The light quarks are present in low energy nuclear physics which is a topic of 
investigation in future chapters of this thesis. 

QCD is a gauge theory based on the non-abelian SU(3) gauge group. The eight 
independent generators of SU(3) give rise to eight massless gluons carrying a color 
charge. Gluons are the strong "force carriers" in QCD. 

The Lagrangian density of QCD is 

Lqcd = IIAF;,F^^'' + q{If - m,)q (1.1) 

where the field tensor F^^ is 

F^u = d,Gl{x)-d^Gl{x)-rigo[Gl{xlGl{x)l (1.2) 

and are the gluon fields. The index a is a color index. The free parameters in the 
QCD Lagrangian density are the gauge coupling constant, go, and the quark masses, 
ruq. 

QCD has been well investigated with perturbation theory in the high energy 
regime. In the low energy limit, QCD should describe nuclear physics and the hadron 
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mass spectrum. Hadron masses depend on the gauge coupling constant like Mhadron ~ 
g-i/so When the QCD coupling constant is large, perturbation theory is not valid and 
a new recipe is needed. The only solution in present day physics is Lattice Quantum 
Chromodynamics (LQCD). Lattice QCD was first introduced by Kenneth Wilson in 
1974 [2]. 

Lattice QCD implements field quantization through path integrals and the dis- 
cretization of space-time onto a four-dimensional Euclidean lattice. The path integrals 
on this space-time lattice allow the lattice gauge theory to be studied numerically with 
Monte Carlo simulations. These simulations share similarities with statistical models 
in Solid State physics. These similarities allow the particle physicist to use similar 
analysis techniques as used in the Solid State models to extract meaningful results 
from the lattice. 

The strangeness contribution to the electric and magnetic properties of the nu- 
cleon has been under investigation experimentally for many years. Lattice calculations 
of the strange quark in the presence of a nucleon are both computationally expensive 
making meaningful results difficult to extract. New computational and numerical 
techniques are needed to determine the nucleon properties. In this dissertation we 
will discuss effective methods using LQCD that will lead to a better understanding 
of the strangeness contribution to the nucleon. 

1.1 Experimental Motivation 

More accurate theoretical predictions of the disconnected strangeness matrix 
elements are needed to compare with experiment. The current experimental mea- 
surement of the low-momentum transfer of the strange nucleon form factors are be- 
ing conducted by groups at HAPPEX [3], A4 [4], and SAMPLE [5]. The most recent 
experimental results published by these groups and the group at Thomas Jefferson 
National Accelerator Facility (JLab) are summarized in Fig 1.1. 

Figure 1.1 is a plot of linear combinations of the electric G's(g^) and magnetic 
Gm(?^) form factors using a parity violating electron- proton scattering process. The 
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Figure 1.1. Experimental results for the simultaneous strange electric and magnetic form 
factors at small four-momentum transfer. 

ellipsed region is the experimental 95% confidence region. The leading lattice result 
marked as [21] is well within this 95% confidence region indicating small positive 
values for the electric and magnetic form factors. Result [21] from the authors of 
reference [1] are from a quenched lattice calculation employing chiral perturbation 
models to extend to the continuum theory. This result can be improved by introducing 
smaller quark masses to the simulation so that a stronger connection can be made 
with chiral models. The agreement with experimental results is strong motivation to 
look deeper into the strange disconnected form factor. 

To make better connection with these experimental results smaller quark masses 
must be used in the lattice calculation. The Wilson QCD action can suffer from 
gauge configurations which produce unphysical results that prohibit the calculation 
of small quark masses. Therefore, theorists must turn to other methods that can 
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avoid these types of damaging configurations. One such method that removes the 
unphysical results and produces more rehable physics is twisted mass QCD (tmQCD). 
The tmQCD action is used in this thesis to improve the strangeness calculation. 

In this thesis, we will discuss the basic lattice techniques that are used in this 
hadron calculation. In chapter two, the basics of lattice gauge theory are reviewed. 
Next, there is a review of twisted mass LQCD and the symmetries that are preserved 
in this formalism. We consider the lattice techniques necessary to extract meaningful 
results in chapter four. 

New work is presented in chapter five. This work discusses new mathematical 
algorithms to efficiently solve linear systems of equations giving quark propagators 
for both the Wilson and twisted mass formalism. In addition to these new methods, 
a perturbative method to calculate the strange quark vacuum expectation values is 
discussed in chapter six. Here, an extension to the existing method in reference [6] 
is employed and an introduction to a twisted mass disconnected perturbative tech- 
nique is given. Finally, the simulation details and numerical results are presented in 
chapter seven. Conclusions of the strangeness calculation and plans for future work 
are summarized in chapter eight. 



CHAPTER TWO 
Lattice Gauge Theory 

Lattice gauge theory is the discretization of the QCD action onto a four- 
dimensional hyper-cubic lattice with a finite lattice spacing. There are, of course, 
an infinite number of ways to define a discrete gluonic and fermionic action on the 
lattice but the simplest method is the Wilson gauge action using the Wilson Dirac 
operator. These methods retain the necessary symmetries that continuum QCD re- 
quires. In this chapter the fundamental concepts of lattice gauge theory are discussed. 
A more complete discussion of lattice gauge theory can be found in many texts and 
journals [7-12]. 

2.1 Lattice Gauge Fields 
The continuum gauge fields are represented by A^, which belong to the gauge 
algebra. The corresponding lattice gauge fields, U^{x)^ belong to the the gauge group 
G. The role of the lattice gauge fields is to move color locally between nearest neighbor 
lattice sites. On any plane of the lattice we define two unit vectors fi and z> that define 
the directional orientation of the gauge links (See Figure 2.1). 



Figure 2.1. A plane in the lattice showing the gauge link structure. 
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Let a be the lattice spacing. If U^{x + a) is the gauge hnk between space-time 
points Xi and Xi + ain the direction, then the gauge field that moves in the opposite 
direction from Xi + ato Xi is the Hermitian conjugate of Uf^{x + a) due to the unitarity 
of gauge fields. 

The continuum and lattice gauge fields are related by 

U^{x) = 6-^"^"^"^^) (2.1) 

where a is the lattice spacing, Qo is the coupling constant and Af^{x) are the continuum 
gauge fields. Gauge fields on the lattice must obey local gauge transformations as 
they do in the continuum theory. To apply a local gauge transformation to a link we 
must specify a gauge transformation at the beginning and end-point of that gauge 
link. Let the local gauge transformation be G{x). The gauge link and fermion fields 
under a local gauge transformation G(a;)are 

U^{x) ^ Gix)U^G''ix + a^), (2.2) 

^(x) ^ G{x)^{x). (2.3) 

With these definitions we are now able to construct gauge invariant operators on the 
lattice. For example, in the pure gauge theory it is now possible to construct a closed 
Wilson loop. A Wilson loop is constructed by taking the trace of four links around a 
closed loop in the - u plane. This operator is independent of starting position and 
is invariant under gauge transformations. The simplest non-trivial Wilson loop is the 
average plaquette. A plaquette is a closed loop, gauge invriant object constructed of 
gaugelinks on the lattice. The average plaquette is an order parameter of the Wilson 
theory. 

According to Wilson, the discrete gauge field action is given by 

where the sum is over all elementary plaqucttcs, U (p). Wilson showed that this action 
is equivalent to the continuum action to leading order in the lattice spacing a. 



2.2 Lattice Fermions 
The Euclidean continuum fermion action for QCD is 



Qcont. 

J p — 



(2.5) 



The four components of If are the usual If = -0^7/^. The 7^ matrices are a 
set of four matrices that satisfy the algebra 



We also define the quantities 



75 = 71727374, 



75 = 75- 



The representation for the 4x4 7 matrices we use is 



a,; 



Ti 



,74 



1 
-1 



(2.6) 



(2.7) 



(2.^ 



75 



v 



-i 

1 



where the index i — 1, 2, 3 and the ctj are the 2x2 Pauh matrices. 

A discrete representation of equation (2.5) is needed for lattice calculations. We 
require that the fermion fields and operators only exist on the lattice sites themselves. 
This is in contrast to the hnks that only exist between lattice points. The lattice 
fermion fields are Grassmann- valued fields that carry flavor, color, and Dirac indices. 
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2.2.1 Naive Fermion Action 

Lattice fermions in Euclidean space are represented by anticommuting Dirac 
spinors, ip{x), that satisfy the relations 

[V',V']+ = [V'^V']+ = [V'^V'1+ = o. (2.9) 

To find a discrete fermion action for these fields, Wilson replaced the covariant deriva- 
tive in the continuum action with a symmetrized difference equation. By using the 
correct choice for gauge links as well, the discrete fermion action remains gauge in- 
variant. To leading order in a, the naive action for the fermion fields is 

S^''^'"^ = mqJ2^i^)^i^) (2-10) 

X 

+ -^Yl '^i^ht^[Ut^i^)'^i^ + )") - Ulix - ii)ip{x - 11)] 

X 

= Y.i^{x)M^^^^^m{y) (2.11) 

X 

where the naive interaction matrix is 

Mi;--[t/] = m,5,y + ^ 5^7.[^.,.<^.,.-M - Ul.,,y5x,y+,]. (2.12) 

A* 

In equation 2.12, ruq is the quark mass and the sum is over Dirac indices. The naive 
fermion action creates huge problems on the lattice. Consider the inverse of the free 
field propagator in momentum space: 

S-^{p) = J]M^^f''^[t/= l]e^^"(^-^). (2.13) 

= mq + -Y'y^sin{p^a). (2.14) 

In the limit as mq — > 0, the inverse propagator creates 2^ zeros in the momentum space 
unit cell. Each of these zeros corresponds to a species of fermion on the lattice. This 
is obviously an unacceptable result. This phenomena is known as fermion doubling 
because there are two species in each direction of the lattice. 
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2.2.2 Corrected Fermion Actions 

There are many possible corrections to the naive fermion action that will remove 
the doubling problem and still remain a "good action" in the continuum limit. Three 
good choices for actions are the Wilson, Kogut-Susskind, and twisted mass fermion 
actions. The advantages and disadvantages of each of these actions will be presented. 

2.2.3 Wilson Fermions 

One approximation to the naive action is the Wilson fermion action. Wilson 
added a second derivative term to the naive fermion action that results in a rescaled 
factor that is related to the bare quark mass by 

1^=^—^ ^. (2.15) 

K is known as the hopping parameter. Equation (2.15) can be solved for the quark 
mass ruq in terms of lattice parameters n and r. The quark mass then is 

mga = 4r (2-16) 

Kc— 1/8 for the non-interaction theory. The same formula holds for the interaction 
case. 

This discrete fermion action, also known as the Wilson action, is written 
= K^i){x)il^{x) (2.18) 

X 

+ ^ X^[^(a;)(7M - r)U^{xyip{x + /x) - ^(a;)(7^„ + r)Ul{x - fi)ij{x - /x)]. 

In the free field limit, when r = 1 the doubling problem is resolved. The matrices 
(7^ — l)Ui^{x) and (7^ + i)Uj^ in (2.18) are the forward and backward quark hopping 
terms, respectively. 

We can rescale the fields in the Wilson action by letting tjj — > \/2K,tjj, giving a 
convenient form of the Wilson action 
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= ^V'(x)V^(x) (2.19) 

X 

+ - r)U^{x)ip{x + ^{x){-f^ + r)Ul{x - ij)ip{x - jj)]. 

It is known that for small quark mass, a; ruq Kc — k. By definition, Kc 
in (2.16) is the value which causes the pion mass to be zero. The calculation of Kc is 
statistical in nature and is determined by the limit ttItt 0. When a zero mode occurs 
at a value of «; < for a given configuration, the quark propagator becomes singular 
in a physical region. These unphysical modes are called "exceptional configurations" , 
and are a large concern for the Wilson action in the quenched approximation (see 
section 2.3). Dealing with this problem is a major focus of this thesis. 

A consequence of the "r" term in the Wilson action is that it breaks chiral 
symmetry at 0{a) in lattice spacing. Consequently, an additive mass renormalization 
is required. The loss in chiral symmetry results in operator mixing and additional 
field renormalizations. 

Even though the Wilson action introduces "exceptional configurations" and 
breaks chiral symmetry, it does preserve a one-to-one correspondence between the 
Dirac and flavor degrees of freedom and the continuum theory. This is a huge advan- 
tage because it allows the interpolating field operators to be constructed in the same 
manner as in the continuum limit. For example, ip{x)ip{x) (scalar) and ip{x)jnip{x) 
(vector) have the same form on the lattice as in the continuum. 

An alternative formalism that is closely related to the Wilson action is the 
twisted mass action. In this formalism an additional term is added to the Wilson 
action that removes the unphysical quark modes. This formalism was proposed by 
Prezzotti and Rossi in 2001 [13]. Twisted mass LQCD is the new frontier for lattice 
calculations and will be discussed in depth in future chapters. 
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2.2-4 Staggered Fermions 

Staggered Fermions reduce the number of fermion species by using one compo- 
nent "staggered" fermion fields rather than the usual four component Dirac spinors 
and by employing a spin diagonalization of the spin components of the fermion 
fields [14-16]. Each of the staggered flavor and spin fields is placed on a corner 
of the lattice. The diagonahzation of the fermion fields removes the 16-fold doubhng 
problem of the naive fermion action. This discretization of the action also preserves 
chiral symmetry when ruq — > 0, because there is no rotation under the subgroup U (1) 
from the single spin index. When chiral symmetry is desired, staggered fermions are 
preferred to Wilson fermions. The exceptional configuration problem is also greatly 
reduced and one can go lower in quark mass in computer simulations. 

The disadvantage of this formahsm is that there is now a 4-fold degeneracy 
for each physical flavor in the continuum limit. The degenerate states are called 
"tastes", to distinguish them from the physical flavors. This degeneracy breaks the 
flavor symmetry at 0{a) on the lattice which makes construction of operators with 
correct quantum numbers difficult. Computationally, staggered fermions save roughly 
a factor of 4 in computer time because they use only a single component Dirac spinor, 
thus saving on storage space as well. 

2.2.5 Lattice Errors 

In any lattice calculation there are statistical and systematic errors. The sta- 
tistical errors are a result of the Monte Carlo stochastic method and fall off like -4^. 
The systematic errors are a result of approximating a spatially and temporelly infi- 
nite problem on a finite lattice. Two well known errors that are a direct result of the 
discretization of the lattice arc the finite volume and finite lattice spacing effects. 

Another source of systematic error occurs when the lattice results are extrap- 
olated to the continuum limit. One must implement a chiral perturbation theory to 
reach the continuum. This extrapolation carries inherent error that appears in the 
final lattice result. 
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2.2.6 Finite Volume Effects 



The volume of the lattice is given by 



Viat — Lx * Ly * Lz * L-, 



(2.20) 



where = arii. Hi is the number of lattice sites in the i direction. If Lj is large, it 
has been shown that the finite volume errors fall off exponentially [17] 



To avoid finite volume effects the length of the lattice must be larger than the particle 
cross-section. A hght hadron cross section is about 2 fm in diameter. Since the 
lattice employs periodic boundary conditions the hadron on the lattice will also have 
reflections of itself in any given periodic direction. When Lj is large enough the 
hadron does not overlap with its reflected image and the volume effect is small. On 
the other hand, if the lattice length is smaller than the hadron diameter and the 
hadron overlaps with it's image, the hadron mass will be large. This produces large 
flnite volume errors. 

2.2.7 Finite Lattice Spacing Effects 

Fields in quantum theories suffer from ffuctuations at all length scales. In per- 
turbation theory, these ffuctuations are responsible for ultraviolet sensitivities and 
infinities in loop diagrams. In light of this, it is hard to understand how we might 
define a discrete approximation to a continuum field that is already randomly fiuc- 
tuating and coarse. Fortunately, only long wavelength objects are physical on the 
lattice. In general, any low momentum, long-wavelength probe is only sensitive to 
space-averaged fields on the order of the probe itself. The averaging of the fields 
suppresses the quantum ffuctuations on the lattice. Consequently the infrared behav- 
ior is not sensitive to a speciffc ultraviolet theory. There are, therefore, an inffnite 
number of ways to construct an ultraviolet theory with the same infrared physics. 



error {rriT^) — e 



—rrinLi 



(2.21) 
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In quantum theory the infrared modes can be affected by the quantum fluctua- 
tions of the ultraviolet mode via the mass and coupling terms. However, if we choose 
an ultraviolet theory that permits us to change the bare coupling and mass terms such 
that the infrared behavior is the same in the continuum limit up to a renormalization 
of 0(a), we can avoid quantum fluctuations [17] . Effectively, the lattice acts as an 
ultraviolet cut-off that restricts the particle modes to low momenta. Ultimately, to 
avoid quantum fluctuations and costly renormalizations in a lattice measurement, the 
lattice spacing a must be smaller than any important scale for the hadron calculation 
under investigation. 

2.2.8 Chiral Extrapolations of Light Quark Masses 

The quark masses, u and d, are too light to simulate in current lattice cal- 
culations because of the exceptional configuration problem and increased statistical 
fluctuations. While new methods are being formulated, the lowest pion mass that can 
be calculated is approximately 500 MeV for the Wilson formahsm. (One can go much 
lower with staggered fermions, but there are interpretational problems.) The current 
method to determine the physical pion mass is to calculate many different pion masses 
and extrapolate to the physical value near 140 MeV. This extrapolation process is 
known as Chiral Pertrubation Theory (xPT). As with any statistical measurement, 
the extrapolated physical m^r has an associated uncertainty. This technique has pro- 
vided reliable results for many lattice calculations, however, the ultimate goal is to 
produce better simulations of the light quark masses so that there is less dependance 
on xPT. 

2.3 Quenched Approximation 
Full QCD calculations are currently unrealistic computationally. A remarkably 
good alternative to full QCD is the Quenched QCD (QQCD). It consists of neglecting 
the determinant of the quark matrix in the lattice gauge field action. Physically, the 
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quenched approximation is equivalent to neglecting the vacuum polarization effects 
of quark loops in lattice calculations. Neglecting these vacuum loops only changes 
the relative weighting of the background for QQCD. 

At short distances the only difference between quenched and full QCD is a small 
change in the QCD coupling constant. This is known as asymptotic freedom. The 
quenched approximation saves factors of 10^ - 10^ in computer time while preserv- 
ing asymptotic freedom, confinement, and the chiral symmetry breaking that QCD 
includes. All of our calculations are performed in the quenched approximation. 



In practice, to generate gauge fields for Lattice QCD Monte Carlo methods are 
employed for the numerical integration of Feynmann path integrals. Monte Carlo 
methods are especially useful in studying physical systems with a large number of 
coupled degrees of freedom in which the inputs have significant uncertainty. 

The QCD path integral is 



where the integration is over gluonic and fermionic fields. The associated QCD action 
with this path integral is 



where M is the fermion matrix. 

As an instructive, simple example [18], consider the path integral of a particle 
moving in a one dimensional well 



2.4 Gauge Field Construction 




(2.22) 




(2.23) 




(2.24) 



where the discrete action is 
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Figure 2.2. Classical particle one-dimensional trajectory (smooth and discrete) from Xi — > 

Xf. 

The corresponding picture of this action is in figure 2.2. 

For large values of A^, the path integral can be determined using a Monte 
Carlo method. A set of possible {x'^s} from i = 1,...,N is a configuration. The 
exponent of the action in the path integral is analogous to the Boltzmann factor in 
statistical mechanics and, thus, is the weight for generating a specific configuration. 
To maximize the efficiency of the method, we wish to generate configurations weighted 
by e""^. This process is known as importance sampling. 

A method that uses importance sampling is the Metropolis method. This 
method begins with an initial configuration and then slightly perturbs each Xi of 
that configuration by a small, random number. This gives a small change in the 
action, AS. After the perturbation, if AS < then the change to the action is ac- 
cepted, otherwise another uniformly distributed random number is generated and the 
procedure is repeated. Each iteration of this method is known as a sweep. To insure 
statistical independence many sweeps occur between accepted configurations. Per- 
forming the Monte Carlo method iterations to obtain independent field configurations 
is called thermalization. 
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A set of configurations is an ensemble. Calculations on the lattice can then be 
performed using the ensemble of the configurations. For the one-dimensional particle 
in a potential well we can calculate, for example, the quantized energy levels of the 
particle can be determined. 



CHAPTER THREE 



Twisted Mass QCD 



As discussed in previous chapters, Wilson fermions are a good solution to 
the fermion doubling problem but introduce zero quark modes which correspond 
to massless quark flavors that produce large, unphysical statistical fluctuations in the 
quenched approximation. A solution was proposed by Prezzotti etal. in 2001 that 
removes the exceptional configurations while retaining the original Wilson symme- 
tries [13]. It is called twisted mass QCD (tmQCD). 



A conceptual problem arises for Wilson fermions in the quenched approxima- 
tion. As we know from field theory, the fermionic determinant contains information 
about the vacuum polarization loops. The quenched approximation neglects the vac- 
uum loops and thus the fermionic determinant. When the determinant is removed, 
exceptional gauge field configurations occur, resulting in large statistical fiuctuations 
leading to a corrupt ensemble average [19]. There have been several regularization 
of the Wilson action schemes proposed to solve this "exceptional problem" [20-22]. 
However, this problem is common to all lattice regularizations using Wilson fermions. 

One solution to the "exceptional problem" is to add a non-standard mass term 
to the Wilson quark action. The lattice Dirac operator is then 



where Dw is the massless Wilson Dirac operator, is the bare quark mass, /x^ is the 
twisted mass parameter, and is the third component of the Pauli matrix acting in 
isospin space. The lattice tmQCD action is then. 



3.1 Introduction to Twisted Mass 



DtmQCD = Dw + mg + iiJ.q-f5r^, 



(3.1) 




(3.2) 



X 
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The tmQCD term in (3.2) generalizes the Wilson fermion action by introducing 
a chiral phase between the mass and Wilson term [23]. As stated above, the twisted 
term protects the tmQCD action from zero quark modes. The protection that the 
twisted mass action offers can be seen explicitly by manipulation of the determinant 

of DtmQCD 

< Det[Dw + mq + inq^y^r^] 

= det[{Dw + mg)\Dw + rUg) + nil (3.3) 

where Det is the determinant in two-flavor space and det is the determinant in one- 
flavor space [24,25]. If the twisted mass term is non-zero, the determinant in 3.3 
can not be zero thus avoiding zero quark modes. Numerical evidence is provided in 
reference [26]. 

The twisted mass parameter couples to terms in flavor space and protects the 
Dirac operator from zero quark modes [13]. Two distinct twisted mass flavors are 
generated from this Dirac operator corresponding to the elements of r^. The twisted 
mass term associated with +1 is the "up" flavor. Likewise, the term associated with 
— 1 is the "down" flavor. To avoid confusion with the up and down quark flavors the 
twisted flavors will be denoted "tmU" and "tmD" for clarity. 

3.2 Classical Continuum Theory 
The continuum twisted mass QCD action is, 

SF[i'{x)ij{x)] = - j d''xij{If + m + i/xg75r3)^. (3.4) 

The axial (75) transformation of the fermion fields is 

which leaves the twisted action invariant [13] and transforms the mass parameters 

m' — mcos^a) + fiqsin{a), (3.6) 
/i' — —msin{a) -\- iiqCos{a). (3.7) 
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where one defines the rotation angle of the transformation by 

tan(a) = ^. (3.8) 
m 

Notice with this definition of the twist angle the standard action is obtained when 

The chiral symmetry of the massless action defines and to be 

= V^T/.TsyV', (3.9) 
= i'lj-^i^. (3.10) 

It is important that the usual symmetries continue to hold in this formalism. 
At non-zero quark mass, the partially conserved vector and axial relations (PCVC 
and PCAC) take the form 



d^Al = 2mP" + (3.11) 
d^V^ = -2/.,e3-^P^ (3.12) 

where the pseudo-scalar and scalar densities are defined to be 

^" = VS75yV','5° = ^V'. (3.13) 

The transformation of the quark and anti-quark to the primed basis results in a 
transformation of the usual Wilson operators. Useful examples of this transformation 
are seen in [13]. The axial and vector currents in the primed basis that utilize fields 
from (3.5) are 

= V^T/.TsyV'' = cos{a)Al + e^^hin{a)V'^ (3.14) 
F;« = V^7^yV'' = cos{a)V^ + sin{a)Al, (3.15) 

for a = 1, 2. When a — ?> these currents have the form 

= VS7/.75y^' = 4 (3-16) 
V'^^i^lJ-^i^' = V^. (3.17) 
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Similarly, the pseudo-scalar and scalar operators in the primed basis are 

P'" = P",(a = l,2) (3.18) 
P'° = cos(a)P^ + ^sm(a)5°, (a = 3) (3.19) 
S"^ = cos(a)S'^ + 2isin(a)P^, (a ^1,2,3). (3.20) 

It is important to notice that in general there is mixing between the axial and vector 
currents as well as the pseudo-scalar and scalar densities. Using the rotated masses 
defined in (3.6) it can be shown that PCAC and PCVC relations take their usual 
form in the primed basis, 

d^A';^ = 2m'P'" (3.21) 
= 0' (3-22) 
with the requirement that the rotation angle is defined as in 3.8. 

3.3 Symmetries of the Bare Theory and Renormalizability 
The massless Wilson Dirac operator in equation (3.2) is 

1 ^ 

Dw = ^ E(^^(^M + V;) - aV;V^). (3.23) 

fj,=0 

The massless Wilson Dirac operator is not invariant under a left multiplication 
of the axial rotation in (3.5) and therefore the Dirac operators are different when 
IJ>q ^ and /iq — 0. This is a welcomed consequence because the twisted mass term 
in the axial rotation protects the action from zero quark modes. If this were not the 
case, the tmQCD theory would still suffer from "exceptional configurations" . 

It has been shown that in tmQCD there is a C/(l) flavor symmetry that leads to 
conservation of fermion number. A vectorial U{1) isospin symmetry also exists which 
is generated by y- 

The twisted mass lattice action is invariant under axis permutations. However, 
reflection symmetries, such as parity, are a good symmetry only in combination with 
a flavor exchange between "tmU" and "tmD" 
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(3.24) 

which is the equivalent to changing the sign of the twisted mass parameter /ig — > —/iq. 
This is a P X ti,2 symmetry of the twisted action. 

Lattice symmetries and power counting prove that the tmQCD model is renor- 
malizable [27]. The P x ti_2 symmetry rules out odd parity, pure gauge terms propor- 
tional to tr[FF] as a —> to contribute to the action [28]. While the couphng constant 
gl and the twisted mass parameter /ig only require a multiplicative renormalization, 
the bare quark mass m needs an additive and a multiplicative renormalization. 

The relationship between the bare and renormalized action parameters are 

gl = Zg{gl,amg,ang,an)gl, (3.25) 
I^R = Zf,{gl,amg,ang;afj.)fj.g, (3.26) 
rriR = Z^{gl,amg,ang]an)mg, (3.27) 

where the Z's are the renormalization factors. The renormalization factors can be 
written in a mass-independent scheme and can be chosen to be independent of arUg 
and a/iq [24]. The mass-independent renormalization parameters are obtained by 
renormalizing in the chiral limit [24, 29]. 

gl = Zg{gl;aij)gl (3.28) 
/^R = Zi^(go;ai^)l^q, (3.29) 
mR = Z^{gl; aii)mg. (3.30) 

Assuming that the massless Dirac operator in (3.23) is of 0(a), then the 0(a) 
improved bare parameters of the action are 

gl ^ glil + bgarrig), (3.31) 
rriq mg + hmam^ 6^0/^5, (3.32) 
jjLg iJ,g{l + bf,amg), (3.33) 
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where nig is the difference between the bare mass and the critical mass, = 
rrio — TTicriticai- The improvement coefficients for the renormahzation bi^,bjn,bfn,bg 
are determined by perturbation theory as well as the PCVC and PC AC relations 
for tmQCD [30]. 

34 TmQCD at Maximal Twist 
Recall that the twist angle defined by the field transformation is defined to be 

tan(a) = ^. (3.34) 
m 

Two interesting choices of the twist angle are a = and a = |. Assignment of a zero 
twist angle returns the standard Wilson lattice action. Choosing a twist angle ct = | 
causes the mass, m, to vanish and is referred to as a maximal twist value. 

As seen in 3.14 a generic rotation by a mixes the axial and vector currents. 
However, when we choose the maximal twist value, there is no mixing but the role of 
the vector and axial currents are exchanged. 

There arc many possible definitions of the maximal twist value. One possibility 
is the Wilson definition of maximal twist. The twist parameter is determined by the 
standard Wilson action when a — 0. The pseudoscalar meson (pion) is calculated 
as a function of the critical mass (hopping parameter k,) and then extrapolated to 
vanishing pion mass. The critical mass parameter is [23] 

arric — 4. (3.35) 

ZKc 

The Wilson definition of maximal twist has been used in previous calculations [31, 32]. 

The tmQCD action expressed in terms of the twisted fields (3.5) has a parity 
violating mass term. This mass term may be removed by a field redefinition where 
the parity violation is now associated with the Wilson term. The resulting action is 
said to be in the physical basis [33]. The parity conservation definition of the twist 
angle is found by enforcing the physical property that there should be no mixing of 
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the charged psuedoscalar and vector current in the physical basis [23, 34, 35], 

J2<V-{x,t)P+{0)>^0, (3.36) 

X 

where the charged pseudoscalar is 

P+{x) = d{x)-f5u{x). (3.37) 

Employing the vector transformation in (3.14) and with the understanding that the 
charged pseudoscalar is invariant under (3.5) we can write the parity definition of 
maximal twist as 

where again the currents with a tilde are constructed in the twisted basis. 

In reference [23] a comparative numerical study between the Wilson and parity 
definitions of maximal twist was performed. Their study showed that there are no 
significant lattice spacing effects on the nucleon or vector meson masses for either 
definition of maximal twist. However, the pion decay constant was found to be 
independent of lattice spacing for the parity maximal twist while the Wilson was not. 

For a fixed value of the twisted mass parameter the parity maximal twist yielded 
smaller pion masses than the Wilson definition. It is desired that the square of the 
pion mass be minimized at maximal twist. The present results imply that the parity 
conserving definition of maximal twist is better for this observable. For this reason, 
the set of {k, /ig) pairs found in [23] will be used in this thesis. 

3.5 Continuum and Chiral Limit 
In tmQCD the lattice cut-off effects resulting from the chiral violating twisted 
mass term may change dramatically as a function of the quark mass. This fact 
is important when chiral symmetry is spontaneously broken. During spontaneous 
symmetry breaking, the chiral phase of the vacuum state in the continuum theory is 
driven by the quark mass term. This is also true in the lattice formalism, therefore 
the continuum limit is taken before the twisted mass /ig ^ [28]. 
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Even with the advancements in computational technologies, lattice techniques 
are not able to compute physical quark masses. Therefore, in the continuum limit, 
a lattice chiral perturbation method is used to reach physical results. Lattice chiral 
perturbation theory (ChPT) is an expansion in powers of the quark mass and the 
lattice spacing parameter that provides estimates of physical observalables in terms 
of a few low energy constants [36] . When ChPT is applied in the tmQCD [37, 38] , 
ChPT involves the renormalized quark mass niR and the rescaled twist angle 

an — tan~^[Ztan{a)], (3.39) 

where Z is a renormalization constants of the operators 'ip'^5Taip' and ■0V' in the mass 
independent scheme described in equation (3.28) in reference [39]. 

0{a) cutoff effects of the pion mass and the pion decay constants are automati- 
cally absent when the twist angle is 90°. However, there are lattice artifacts of O(^) 
that remain from the chiral Lagrangian density in the pion mass [40]. 



CHAPTER FOUR 
Lattice Techniques 

In this chapter a brief review of lattice strategies to extract information from 
lattice calculations is presented. The purpose of this chapter will be to present a 
review of two-point Green function source techniques and correlation functions. We 
will also discuss the strange matrix elements of the nucleon. 

4.1 Grassmann Integration 
Grassmann integration is a useful tool to evaluate fermionic integrals in two and 
three point functions. A brief summary of the properties for Grassmann variables is 
presented here. Let the Grassmann variables and it's conjugate by ( and (*. If these 
are to be Grassmann variables they must obey the anti-commutation relations 

[ao]+ = [c;,o]+ = [c;,c;]+ = o. (4.1) 

Integration over Grassmann variables can be defined as 

dC = f dC = 0, 



1 dCC = / dCC = 1- (4.2) 
Prom equation 4.2 we can deduce the property 

/ U^dCdCmexpi- J2 C:M.,Q] = det{M). (4.3) 

This integral differs from the corresponding integral over commuting variables by 
resulting in the det{M) rather than det{M)~^. 

Suppose now that the Grassmann variables represent a quark field. Then, 
for example, using Wick contractions between quark and anti-quark fields then the 
integral in (4.4) results in a quark propagator 

dCdCU^e-'^^^ = det{M)S^p. (4.4) 
25 
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We set det{M) = 1 for these types of integrals in the quenched approximation. A sim- 
ilar expression can be determined for tmQCD. In chapter 3, the field transformations 
at maximal twist was expressed as 

tptm = ;^(l±«75)V','0tm = ^V'(l±«75), (4-5) 

where the + and — represents "tmU" and "tmD" , respectively. 

We are interested in how Grassmann integration behaves using twisted fields. 
Our example from equation (4.4) using maximally twisted fields can be expressed as 

j dCdCil ± ^75)CaC/3(l ± ^75)e-^'""'^ = det{M'){l ± Z75)^a/3(1 ± ^75)• (4.6) 

with M' = (1 ± i75)M(l ± ^75) and where the propagator, Sap, is in the physical 
basis. Again, let det{M') — 1. This instructive, simple example shows how to create 
quark propagators in the twisted basis and return to the physical basis by twisting 
the ends of the propagator. This strategy was employed to determine hadron masses 
in reference [41]. 

4-2 Green Function Methods for Proton/ Neutron 
In this section we will review the proton two and three point function method 
presented in reference [42] as well as the twisted mass representation. The twisted 
interpolation fields used for the proton two point function are 

Xa{x)tm = 6'^''^Vi"^«(x)i^4")'(x)i^(C)^^Vf^(2^)tm, (4.7) 



tm- 



The Greek and Latin indices represent Dirac and color indices, respectively, in equa- 
tion (4.7). The interpolation fields for the neutron are given by a m — > d field exchange. 

The proton two point function for forward time {t > 0) can be written in terms 
of the interpolation fields as follows: 
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Gpp{t;p,r') = ^e-'^^-X,„<mc|r(xa(x)tmXa'(0)tm)|^ac> (4.8) 

X 

= ^e-^-r:,'ae'^'^(-e"'''^')(C^)/3.(^)y/3' (4.9) 

X 

The 4 X 4 r' matrix determines which correlation function is to be evaluated and is 
generic until specified. A similar function can be written for the neutron using the 
correct interpolation fields; however, we will focus on the proton here for clarity. 

We have defined the charge conjugation matrix C = 72 and C = C75, which 
satisfies the relation C'jf^C'^ = 7*. A general transformation can be constructed for 
a general matrix Q such that Q = {CQC'^Y . 

In Euclidean space the integration formula for the time ordered N-point function 
is defined to be 

< vac\T{xl:^{-UA)i'(}{-itB)...)\vac >= J (it/(iC(iCe"^°"^^'^"^lCa(^A)C/3(tB) 

(4.10) 

where Sq and Sf[C, C] are the Euclidean gluonic and fermionic actions respectively. 

Using Grassmann integration, we may write the proton two point function in 
the physical basis as 

Gpp = Yl e-'^^-'e^''e^'''^{tr[r' + !^^^ ^(")°°'(a;, 0) + !?'^^ 

X V V 

X (1 -;^^) ^(^)'>'>-(^,o) (^ -;^^) +;^^) ^(")--(x,o) (^ +;^^) ] 

+ ^rriilM^Haa'^ o) (l+;75) j^^^ (l-^75) ^(.)..( 0)^1^ 

X li±S^5H-'(x,0)il±S^]), (4.11) 
V 2 V 2 

where a configuration average is understood and the trace is only over Dirac indices. 
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Using the property of traces we can rearrange the multiphcations such that 

X V V 

+ i^[Ii±!75)r'ii±!2:^5H-'(x, 0)]tr[^('^)'"''(x, 0)5(")^^'(x, 0)]). (4.12) 
V 2 V 2 

If we define a new gamma matrix, F*"' = |(1 + i75)r'(l + ^75) it is possible to write 
the proton two point function as 

Gpp = e-'P-'^e^^'e"'"''''' {trp"" S^""^""' {x, Q)S^'^'^^^' {x, 0)5(")"'='(a;, 0)] 

X 

+ tr[r*"'5(")""'(a;,0)]tr[^('^)^^'(a;,0)5(")'="'(a;,0)]). (4.13) 
The proton two point function presented in the Wilson formahsm is [43] 

Gpp = Y e-'^'^e'^^e"'^'"' {tr[TS^''^''''\x, Q)S^'^^^^' {x, 0)5(")'^"'(a;, 0)] 

X 

+ ir[r5^")""'(x,0)]tr[^('^)^'''(x,0)5^")'^'^'(x,0)]. (4.14) 

The form of the two point function is the same in equations (4.13) and (4.14) if 
ptw _ p This discussion shows that the same techniques can be employed as in the 
original Wilson case with the exchange of F — > F'. 

As suggested by [41], in practice the ends of the propagator are twisted upon 
creation of the quark propagators so that calculations can be done in the usual way 
in the physical basis. Since the usual hadronic two-point functions may be used, the 
rest of this chapter will assume we are doing the calculating in the physical basis. 

4-3 Correlation Functions 
Properties of correlation functions are a fundamental concept for analysis of 
hadron structure [42-46]. A review of correlation functions is given in this section. 
In the large time limit {t » 1), the proton two point function is 

Gpp{t;p,r) N^Y^^'^^'^cc'a < vac\xa{0)\p, s >< s,p\xa'{0)\vac >, (4.15) 
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where A^^ is the number of spatial lattice points. In the two-point function we have 
used the fermionic lattice completeness relation 

'^^\n,p, s >< n,p, s\ — I. (4.16) 

n,p,s 

The corresponding continuum completeness relation is 

n,s ^ ' 

Thus, the correspondence between lattice and continuum states is 

\n,P, S >lattice-^ {:^y^^\n,p, S >cont. (4-18) 

where the volume of lattice sites isV = N^a^. With this relation and the continuum 
field relation, t/jiat ^^^^'^^^^'^cont it is possible to determine the matrix elements of 
the interpolation fields in (4.15). These are 

<Vac\Xoc{^)\p,S>iat {2Kf/^ ^l^^^'^ ^ '^Qc|Xa(0)|P,g >cont (4.19) 

<P,s\xa{^)\vac>iat J2^^ijf^f^ <P^^\Xa{^)\vac>cont ■ (4.20) 

(4.21) 

The lattice matrix elements are related to the continuum free spinors Ua{p-i s) and 

UaiP:S) by 

< vac\xa{0)\p,s >iat = Aua{p,s), (4.22) 

< p,s\xa'{0)\vac>iat = A*Ua'{p,s), (4.23) 

where A is a complex scalar in general. Now we are prepared to determine the large 
time limit of the proton two-point function as a function of the momentum and F [42] . 

where the usual relation for free spinor fields has been employed, 

u{p, s)u{p, s) = . (4.25) 
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A similar argument is proposed for the proton three point function. The three point 
function is constructed with a current insertion between the interpolation fields in 
equation 4.8. The large time limit of the three point function is then 

Gpj,,(t2,ti;p,P,r) ^ -^iv2^e-^-(*^-*i)e-^.'*^ X (4.26) 

s,s' 

^a,a' < Vac\Xa{0)tm\p, S > X 

< p, s\Jf,{0)\p', s' Xp', s'\xa'{0)tm\vac > 

where t2 is a time after the current insertion and ti is a time before. Pictorially, the 
two and three-point functions are seen in Figure 4.1. t is the final time index and t' 
is the time step at which the current is inserted in this picture. 



<:t,t-.t:i> — 



Figure 4.1. Two and Three Point correlators. The solid lines represent quark propagators 
and the shaded box is a current insertion. 



The lattice, continuum relation for the current expectation value above is 

< P, s\JM\j/, S' >lat^ ^{^^f^ < p, s\J^iOW, S' >,onU (4-27) 

where the continuum state is 

< p, s\J^{0)\p', s' >cont= iuip, s)(7^Fi - a^y^F2)u{p, s). (4.28) 
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Fi and F2 are real functions and a/^i, — ^[7/^,7;/] 
Given, 



and we choose the zero momentum (p = 0) charge density (/x = 0) as the current, the 
proton three-point function becomes 

G,j,,{h,t,;0,-q,r') - 5e-(*-*^)e-^*n^^)(i^i - (^^2). (4.29) 

2 

Here we identify G(.{q^) = (Fi — ^^^^2 -^2) as the electric form factor of the nucleon. 
Similarly, with 



o-fc 




the zero momentum, space-like {^l — i) three point function becomes 

Gpj,p(t2, h- 0, -g, ^ (*^-*i)e-^*ie,-fe^gK^i + ^2). (4.30) 

For this choice of T we find the magnetic form factor Gm{(f) = {Fi + F2). 

4.4 Strange Matrix Elements 

Once the two and three-point functions are calculated methods are employed 
to extract the electric, magnetic, and strange matrix elements from the correlators. 

A common technique is to create a ratio of the correlators and then sum over 
the time insertion index. The ratio itself is 



Rx{t,t',q) = 



(4.31) 



G'(2)(i,0)G(2)(f,5-) 

where the index X = {E, M, S} are the electric, magnetic, and scalar ratios re- 
spectively. The indices t and t' are the sink (final time) and current insertion time 
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values [47]. The three point function in equation (4.31) is constructed from the cor- 
relation of the two point function with the loop data when a disconnected part is 
evaluated. 

Define the Fourier transform of the self contracted disconnected lattice current, 

J{x, t), to be 

J'(g-',i) = j;e-^-J(f,i). (4.32) 
The disconnected three-point function can then be written generically as [48] 

G^^\t,t',q) = <G^^\t,0)J'{q,t')> (4.33) 
- <G^''\t,0)><J'{q,t')>. 

Strange matrix elements are extracted from equation (4.31). The extracted 
matrix elements are related to the form factors by 

M(.,M,.) = {Gs,'-f^,GE]. (4.34) 

For the magnetic case, k are indices over the spatial directions. All other indices 
for the magnetic form factor are suppressed for simplicity. 

There are many ways to extract the matrix elements from the form factors. One 
way to acquire the matrix element is to sum over the contributions of the inserted 
strange quark currents [49] 

t 

^ Rx{t, t', q) constant + tMx{t, q). (4.35) 
t'=i 

A disadvantage of this method is that it depends on a linear fit of the data, which may 
only be accurate in a specific temporal region [47] . An alternative method employeed 
by reference [50] is 

t fixed 

Rx {t, t',q) ^ constant + tMx (t,q), (4.36) 

t'=i 

where t fixed > t. 
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In both of the previous methods a hnear temporal fit of the data is require to 
measure the matrix element. In practice, the fit is restricted to a limited set of time 
slices. To remove this linear dependance, a differential method can be employed [47, 
48] . Using the form 

t+i 

J2[Rx{t, t', q) -Rx{t- 1, t\ q)] ^ Mx{t, q). (4.37) 
t'=i 

The resulting matrix element is constant over a larger range of time slices and is not 
subject to a linear fit of the data. This method was employed in the high statistics 
study of these matrix elements conducted in reference [47]. 



CHAPTER FIVE 



Linear Equations Solution Techniques 



For either the Wilson or Twisted Mass approach to LQCD, we are faced with 
solving large, sparse systems of linear equations to determine the respective quark 
propagators. This chapter focuses on improving iterative methods for solving these 
systems of linear equations, which often involve multiple right-hand sides and multiple 
shifts. New Krylov iterative methods to solvie these systems of equations will be 
presented in this chapter. 



5.1.1 Eigenvalue Projections 

There are two general types of projection methods used to evaluate eigenvalue 
equations. These two are oblique and orthogonal projection methods. In this thesis, 
we consider only orthogonal projections. Orthogonal projection methods approximate 
an eigenvector 2; by a vector z. 

Let M be an n X n complex matrix and X be an m — dimensional subspace 
of the space C". Our goal is to determine the eigenvalues. A, and eigenvectors, 2;, of 
the eigenvalue equation 



where z belongs to and A belongs to C. 

To determine the projection operator we must find the appropriate eigenpair 
(A, z) for equation (5.1), with A in C and z in i^, such that the Galerkin condition is 
satisfied. The Galerkin condition is the requirement that the vector Mz — XzinK is 
orthogonal to all other vectors v E K, 



5.1 Projection Methods 



Mz = \z, 



(5.1) 



Mz-Xzl. K, 



(5.2) 



34 



35 



which can be written as 

(Mz - ~Xz, v) =0, Vt; e K. (5.3) 

When this condition is true, the approximate eigenvector z is completely contained 
in K and therefore is exact. 

Assume that an orthonormal basis {vi,V2, ■■■,Vm} of K exists and that the 
matrix V is constructed with the vectors Vi,V2, ■■■,Vm as columns. 

In this chapter, (a, b) denotes an inner product between two vectors a and b. 

Let 

z = Vy, (5.4) 

so that equation (5.3) becomes 

{MVy - ~XVy, vj) = 0, j = 1, m. (5.5) 

If we identify the matrix Bm — V^MV, y and A must satisfy 

BmV = Ay. (5.6) 

This provides a numerical method to determine approximate eigenvalues and eigen- 
vectors of M using the Galerkin condition in equation (5.3). This is known as the 
Rayleigh-Ritz procedure and can be summarized in Table (5.1). 

It is possible to reformulate orthogonal projections in an operator language. 
Consider again the Galerkin condition in (5.3). Define the projection operator Pk — 
V'^V . The Galerkin condition becomes 

Pk{Mz -Xz) ^0,XeC,z e K. (5.7) 

Since the operation of the projection operator on the approximate eigenvector z is 
invariant, the operation of Pk on equation (5.3) can be viewed as a linear transforma- 
tion from K to K [51]. Another way to write the operator expression of the Galerkin 
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Table 5.1. Rayleigh-Ritz Procedure 

1. Compute an orthonormal basis {t'i}i=i,...,m of the subspace K. 
Let V — [^1,^2, ^m] whose columns span K. 

2. Compute = V^MV] 

3. Compute the eigenvalues of Bm and select the k desired 
A, i = 1, 2, j where k < m 

4. Compute the eigenvectors y^, i = 1,2, ...,k, of B^ associated 
with X, i = 1,2, k 

and the corresponding approximate eigenvectors of M, 
Zi = Vyi,i = 1,2, 



condition is 

PkMPkz =~Xz,~XeC,zeC'' (5.8) 

which exphcitly shows the hnear operator A^, = PkAPk for the whole space C". If 
we are restricted to an orthogonal space K, this is the matrix B^n- Equation (5.7) is 
known as the Galerkin approximate problem. 

A useful property for estimating the convergence of projection methods for 
eigenvalue equations is the distance || (/ — Pk)z II2 of the exact eigenvector z from 
the subspace K. For this distance we have the inequahty [51] 

\\z-z\\2>\\{I-Pk)z\\2: (5.9) 

such that a good approximation of the eigenvector z from K results when || (/ — 
Pk)z II 2 is small. 

5.1.2 Harmonic Rayleigh-Ritz Procedure 

While Rayleigh-Ritz values do a good job of determining approximate eigen- 
values (Ritz Values) on the exterior of the eigenvalue spectrum, problems can occur 
when interior Ritz values are calculated. When a Ritz value is on the exterior of the 
spectrum, the associated Ritz vector usually has some significance. In contrast, the 
Ritz vector in the interior may be a combination of many eigenvectors in the subspace 
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giving an interior Ritz value with little meaning [52]. These are known as Spurious 
Ritz Values (SRV). Spurious Ritz values can have adverse effects on existing Ritz 
values of significance. When a SRV is near a "good Ritz value" the corresponding 
eigenvectors blend together. In this situation, it is necessary to determine the residual 
norm to distinguish which of the Ritz values is of significance. 

A solution to eliminate the SRV problem is to convert interior Ritz values to 
exterior Ritz values. A modified Rayleigh-Ritz procedure called the 'Interior' or 
'Harmonic' Rayleigh-Ritz procedure is presented [52-54]. The Harmonic Rayleigh- 
Ritz procedure presents a solution to the SRV problem by shifting the interior values 
to the exterior of the eigenvalue spectrum. 

Consider the eigenvalue problem 

Mz = Xz. (5.10) 

Let X be a j-dimensional subspace of C". It is from this subspace that we wish to 
extract the approximate eigenvectors. To extract an interior eigenvalue the shifted 
matrix (M — al)"^ should be used in the Rayleigh-Ritz method. This matrix shifts 
the eigenvalues to the exterior of the spectrum for this operator. The analysis of the 
procedure will make use of this operator, but in practice this shifted, inverted matrix 
is never calculated. Creating this matrix is impractical because of the additional 
computational cost of finding solutions of linear equations. 

Applying the generalized Rayleigh-Ritz procedure to the shifted interior prob- 
lem, we find 

gt(M - aiy'Qd = -^—Q^Qd (5.11) 

U — (TI 

where {9, Qd) is the approximate eigenpair of the matrix M. The matrix Q should 
span the columns of the subspace K. Instead, to avoid having to calculate the in- 
verted, shifted matrix, let Q = {M — a I) P. Equation (5.11) becomes 

P^{M - aiyPd = —^P\M - aI)\M - aI)Pd. (5.12) 
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Solving this generalized shifted and inverted Rayleigh-Ritz equation yields the 
eigenpair (^z^, Qd). This is the corresponding eigenpair for the matrix M. However, 
since we are trying to extract an interior eigenvalue with the shifted, inverted matrix 
(M — al)~^, a better choice for the approximate eigenpair of M is {p,Pd) where 
p is the Rayleigh quotient with respect to M. Pd is a better approximation for 
the interior eigenvector than Qd since we have shifted the problem. Likewise, the 
Rayleigh quotient p is a better approximate eigenvalue of M than 9. 

This analysis has led us to expect that if z is approximately in K, then the 
harmonic Rayleigh-Ritz method will produce a good approximation to z and an as- 
sociated eigenvalue near a. If we let the approximate eigenvalue of the shifted system 
he 9 = a + S, then we may write the harmonic Rayleigh-Ritz equation as 

P\M - aiyPd = \p\M - aI)HM - aI)Pd, (5.13) 



By multiplying by the vector d* and determining the two-norm we find that equation 
(5.13) yields 

\\{M - aI)Qd\\l < \S\\\{M - aI)Qd\\2 (5.14) 
\\{M - aI)Qd\\2 < \S\. (5.15) 

Therefore, if the harmonic Ritz value is within S of the shift a, the residual norm 
must be bounded hy \S\ [55]. For a harmonic Ritz value close to a and in the limit 
5 — > 0, the harmonic Ritz vector cannot be spurious. 

5.2 Projections for Linear Equations 
Projection methods are useful for solving systems of linear equations as well 
as eigenvalue problems. Most practical iterative methods for solving a large system 
of equations employ a projection process at some stage of the algorithm. A few 
good projection techniques that are used are the Galerkin, MinRes, and Left-Right 
projections. 
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5.2.1 General Projection Method for Linear Equations 
Consider the linear system of equations 

M{x-Xo) = ro, (5.16) 

where the nxn matrix M is a complex. Projection techniques are designed to extract 
an approximate solution of the set of linear equations from a subspace of C". Let K be 
the m-dimensional search subspace. There must be m constraint equations to extract 
a solution from the subspace K. The usual way to determine the m constraints is to 
enforce m independent orthogonality conditions. Specifically, we require the residual 
vector r — b — Mx to be orthogonal to m hnearly independent vectors. This set of 
m linearly independent vectors defines another subspace L which is referred to as 
the constraint subspace or left subspace [56]. This general structure is known as the 
Petrov-Galerkin conditions. 

Let the column vectors of the matrix Vnxm = [vi,V2, ■■■,Vm] form a basis for 
K. Likewise, let the columns of Wnxm = [wi, W2, Wm] form a basis for L. If the 
approximate solution vector extracted from the search space is 

x^Xo + Vy, (5.17) 

where Xo is the initial guess, then the orthogonality condition requires that the system 
of equations for the solution vector y must be 

W^fMYy = lyVo. (5.18) 

To is the residual vector associated with the initial solution vector Xo- If the assump- 
tion is made that the mxm matrix W^MV is non-singular, then the approximate 
solution vector is 

x^Xo + V{W^MV)-^W^ro. (5.19) 

The procedure just described is known as the prototype projection method and is 
summarized in Table (5.2). The approximate solution vector, x, extracted from the 
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Table 5.2. ALGORITHM:: Prototype Projection Method 

1. Until convergence, Do 

2. Select a pair of subspaces k and L 

3. Choose bases V — [t"!, t'2, fm] ^-nd W — [wi, ^2, w^] for k and L 

4. r :^ b- Mx 

5. y := {W^MVy^W^r 

6. x := a; + 1^7/ 

7. Enddo. 



Table 5.3. Existence Criteria of the solution x. 

1. M is positive definite and the left subspace L — k or 

2. M is non-singular and L = Mk. 



prototype projection method is only vahd if the matrix B — W^MV is non-singular. 

The matrix B can be singular even when the matrix M is non-singular. If either 
of the following conditions in Table (5.3) hold, then B is non-singular for any bases 
V and W oi K and L and the prototype projection method solution exist [56]. The 
conditions that need to be satisfied are in the Table (5.3). 

Specific projection methods are determined by choosing specific vectors that 
form a basis for the search and left subspaces K and L, respectively. Two common 
projection methods are the Minimal Residual and Galerkin Projection methods. 

5.2.2 Specific Projections for Linear Equations 

The Minimum Residual projection method (MinRes) is created with a specific 
choice for the spaces K and L. For a MinRes projection we choose the left subspace to 
be L = MK. The basis vectors for the left subspace are then W — MV. Therefore, 
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the MinRes projection can be written as 

(MVyMVy = (MVYro (5.20) 
y = {{MVyMV)-\MVyro. (5.21) 

The approximate solution vector is constructed out of the search space as before 
X ^ Xo + Vy. 

The Galerkin projection can be constructed with the choice for the left subspace 

to he L = K . The basis vectors of the left space are W = V. The projected system 
of equations that we wish to solve now is 

V^MVy = {Vyro (5.22) 
y = (0My)-VVo, (5.23) 

where the solution vector is constructed in the same manner as with the MinRes 
projection. Projection methods are incredibly useful in that they project large prob- 
lems of dimension-n into smaller, more manageable problems of dimension-m. This 
property is valuable for many iterative methods discussed in this thesis. 

5.3 Orthogonal Matrices 
In many algorithms an orthogonal basis for the search subspace is needed to find 
a solution to a system of linear equations. A few common methods to create the basis 
vectors of the search space are standard Gram-Schmidt, Householder reflectors, and 
Fast Givens Rotations. We will discuss the numerical advantages and disadvantages 
of these orthonormal rotations in this section. 

5.3.1 Gram- Schmidt 

The set of vectors G — {gi,g2, gr) is said to be an orthogonal set if the inner 
product of all the elements of G are zero when i j. This same set of vectors is said 
to be orthonormal if in addition || gi \\2= 1, Vi. A vector that is orthogonal to all the 
vectors in the set G is said to be the orthogonal complement of G and denoted by 
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Table 5.4. ALGORITHM :: Gram-Schmidt Procedure 



1. 


Compute rii =|| Xi \\2- If rn = 


= Stop, else compute qi = Xi/ru 


2. 


Forj — 2, .., r, £)o 




3. 


Compute = {xj, qi) for i — 


l,2,...,j-l 


4. 


q — Xj i^i^irijqi 




5. 






6. 


If Tjj = then Stop, else qj = 




7. 


Enddo 





G±. A unique vector Xi can be written as the sum of vectors from G and G±. The 
Gram-Schmidt process takes any vector Xj {xj e G±) and orthogonalizes that vector 
with respect to all previous vectors Xi (x, e G) to form an orthonormal set of bases 
vectors. The Gram-Schmidt algorithm is given in Table 5.4. 

Here q is an orthogonal normalization measure in the context of the convergence 
of the algorithm. 

Notice in steps 4 and 5 of the Gram-Schmidt algorithm that the vectors q and 
Tjj are generated with a QR decomposition. A QR decomposition exists whenever 
the column vectors form a linearly independent set. The Gram-Schmidt algorithm is 
a common orthogonalization method, but is known not to be as numerically stable 
as other algorithms. 

5.3.2 Householder Matrices 

An alternative approach to the Gram-Schmidt procedure is to use the House- 
holder algorithm. This technique uses Householder reflectors to build an orthogonal 
matrix. A reflector is a matrix of the form 

g = / - 2ww^ (5.24) 
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where w is a normalized work vector. A reflection matrix that leaves the flrst k — 1 
columns unchanged while zeroing the k^^ column is of the form 



These reflectors geometrically represent the reflection of a vector Xi into some plane. 
To obtain an orthogonal set of vectors using Householder reflectors we construct 



where — {Qm-i---Qi)^ and R is an upper triangular matrix generated from m — 1 
Householder transformations onto X. Householder reflectors have the advantage of 
being more stable than standard Gram-Schmidt but have an additional overhead 
expense due to the multiplication of the work vector w on to itself to form the House- 
holder reflectors. For large matrices, the additional cost of the Householder matrices 
can make the overhead of this algorithm large. 

5.3.3 Givens Rotations 

A fast method that can be invoked to determine orthogonal matrices is the fast 
Givens Rotations [57]. In contrast to Householder Reflectors that eliminate all the 
elements but the first in a given vector, a Givens Rotation eliminates each element 
individually. In a parallel computing environment (such as MPI), both fast Givens 
Rotations and the Householder algorithm can have a significant speed advantage 
relative to Gram-Schmidt. The Householder Method requires 0{nlog{n)) steps and 
{n — 1) square roots using n{n — l) processes while the fast Givens Rotations require 
0{n) steps to create an orthogonal matrix [58]. An example of a Givens Rotation 
matrix is 



Qk = I - 2wkW,^ . 



(5.25) 



(5.26) 
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V 



.. c .. s .. 



.. -s .. c .. 



.. .. .. 1 



where c and s are the Givens cosine and sine, respectively. These trigonometric func- 
tions can be determined exphcitly. For example, to annihilate the bottom element of 
a 2x1 vector we have 



c s 



—s c 



J 





which gives the constraints on the cosine and sine. The constraints are sa + cb = 
and + — 1, which result in the following algebraic form of the Givens cosine and 
sine: 

c = a/^a2 + b^,s = -b/^a^ + b^. (5.27) 



The factorization is then determined by 



Q — GiG2---Gg, 



(5.28) 



where there are g = {2'm+n+l) /2 Givens matrices for a generic m x n matrix M. This 
method to determine orthogonal matrices is preferred when solving large systems of 
equations due to the reduction in overhead of the calculation in comparison with the 
aforementioned algorithms. This method is stable while producing reliable results. 
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5-4 Krylov Subspace Methods 
A general projection method extracts an approximate solution vector from 
the system of equations 

Mx^b (5.29) 

by employing the Petrov-Galerkin condition that requires the residual vector, r = 
b — Mx, to be orthogonal to the left space L. A Krylov method is a method in which 

the test space is a Krylov subspace 

Kra{M, To) = span [vo, Mvo, MV„, M'^-V^} , (5.30) 

where x^ is the initial guess and = h — Mxq. This condition is true for all Krylov 
methods. Krylov methods differ in their choices of the left space and by how the 
problem is preconditioned. It is clear that the approximate solution vectors extracted 
from the Krylov subspace is of the form 

M^^b ^ Xm 

= Xo + qm-i{M)ro, (5.31) 

where g^_i(M) is a polynomial in M generated by Km{M,ro). The choice of the 
left space, which is generated by the constraints used to build these approximations, 
will have an important effect on the particular iterative method. Two examples of 
choices are for the MinRes projection in which Lm — MK^ and the Galerkin 
projection with = K^. 

5.5 Arnoldi Method 
Arnoldi's method is an orthogonal projection method onto a Krylov subspace 
of dim{m) for general non-Hermitian matrices. The Arnoldi procedure can be used 
both to compute eigenvalues and to solve systems of linear equations. The Arnoldi 
procedure to build an orthogonal basis is listed in Table 5.32. At any step in the 
algorithm the previous Arnoldi vector, Vj, is multiplied by the matrix M to form 
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Table 5.5. Arnoldi Algorithm 



1. Choose a vector vi such that || vi ||2 = 1 

2. For j = 1,2, ..,m, Do 

3. Compute hij = {Mvy, vi) for i = 1,2, 

4. Compute Wij := Mvj — T,l_^hijVi 

5. hj+ij =\\wj\\2 

6. Ifhj + 1, j = then stop 

7. Vj+i = Wj/hj+ij 

8. Enddo 



This vector is orthonormalized against all previous Vi vectors with a standard 
Gram-Schmidt procedure. The set of vectors, form an orthonormal basis of 
the Krylov subspace. Let Vm be the n x m matrix whose columns are {vi...Vm}- 
Let Hm be the m x m upper-Hessenberg matrix formed by the hij values from the 
algorithm. Then the Arnoldi iteration gives the recurrence [56] 

MVm = VmSjn + ^m+l,m'^m+lC^ = Vm+lH,^, (5.32) 

giving, 

V^MVm = H^. (5.33) 

Pictorially we can see the action of M on the basis vectors Vm in Figure (5.1). 
We first consider how the Arnoldi recurrence can be used for eigenvalue computa- 




Figure 5.1. The action of M on Vm gives VmHm plus a vector 
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tions. Essentially the Arnoldi algorithm combines use of a Krylov subspace with 
the Rayleigh-Ritz projection. Arnoldi concluded that the eigenvalues of a Hessen- 
berg matrix smaller than the dimension of the original matrix can provide accurate 
approximations to some eigenvalues of the original n x n matrix [56]. Once these 
approximate eigenvalues are known, an approximate solution to the original problem 
can be determined. 

As a result of the projection onto Km we gain the approximate eigenvalues A^-"^^ 
of the Hessenberg matrix Hm [51]. The approximate eigenvector associated with the 
the eigenvalue Xf^^ is defined to be 

dt^^VmvtK (5.34) 

Using the eigenvalue equation, the small i*'* eigenvalue problem is then 

V'^MVdi = Oidi, (5.35) 

where 9i is the i*'* approximate eigenvalue. The associated Rayleigh-Ritz approximate 
eigenvector is i/i — Vdi. The eigenvectors and eigenvalues form Rayleigh-Ritz pairs 
{Qi-iUi) where t/I™"^ is the associated eigenvector of length m. 

For a moderately sized Krylov subspace, a few of the approximate eigenvalues 
are usually good approximations to the true eigenvalues of the original matrix M. 
As the dimension of the Krylov subspace increases, the quality of these approximate 
eigenvalues improves until all of the desired eigenvalues of M are found. Obviously 
it is not practical to have a large Krylov subspace due to storage and computational 
cost. However, with a reasonable dimensioned subspace, the Ritz eigenvalues can 
play an important role in deflated Krylov methods. It is important to be able to 
cost-effectively estimate the residual norm during Krylov method iterations. A cheap 
way to determine the residual norm makes use of the expression [51], 

(M - A^)/)«l") = hm^,,melyt^vm^,. (5.36) 
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The two norm of equation (5.36) is 

II (M - Xt-h)ut^ ||,= h^^,,^ I e^yt^ \ . (5.37) 

So, the residual norm is equal to the last component of the eigenvector y^^^ multiplied 
by hm+i,m [51]. 

When multiple shifts are desired with Krylov methods, the Arnoldi iteration 
can be modified to handle these shifted systems of equations [52]. 

The new shifted operator (M — al) gives the eigenvalue equation 

V^{M - alf{M - aI)Vd ={6- aI)V^{M - alfVd, (5.38) 

where 6i is a Harmonic Ritz value. The harmonic Rayleigh-Ritz pairs are {9i, iji) where 
we have used the relation y^ — Vdi. When the harmonic Rayleigh-Ritz procedure is 
applied to the Arnoldi iteration we have the eigenvalue equation [52, 59-62] 

(H^ + hl^,^^fel)d^ed, (5.39) 

where / = {Hm — al)~^ejn- The problem has now been altered from finding eigen- 
values and eigenvectors of Hm to finding the eigenpairs of equation ( 5.39). 

5.6 Arnoldi Methods for Linear Equations 
We next consider using the Arnoldi recurrence for solving linear equations. 
These are ways of applying the projection techniques from section 5.2 to a Krylov 
subspace. The choice of the left subspace determines the iterative technique. In the 
next section we will introduce popular methods that are widely used for different 
choices of the left subspace L^. 

5.6.1 Full Orthogonalization Method 

We consider an orthogonal projection method for a system of equations Mx — ro 
which uses the left space = Km = Km{M,ro), with Km defined in equation 
(5.30). This method finds an approximate solution vector Xm from the affine subspace 
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Xo + Km{M, To) by imposing the Petrov-Galerkin condition 

b - Mxm ± Km. (5.40) 

If the first basis vector of the Krylov subspace in Arnoldi's method is vi — 

ro/\\ro\\2, then 

V^MVm = Hm (5.41) 
holds with (3 = ||ro||2. If we then employ (5.32) we may write 

VZro^V^{Pv^)^Pe^. (5.42) 

This results in the approximate solution vector 

Vm = H-\Pe,), (5.43) 

Xm = Xo + VmVm- (5.44) 

The Arnoldi method for linear equations with a Galerkin projection is referred to 
as the Full Orthogonalization Method (FOM) [56]. The FOM algorithm is described 
in Table 5.6. 

5.7 GMRES Methods 

5.7.1 Standard GMRES 

The General Minimal Residual method (GMRES) is the MinRes projection 
applied to a Krylov subspace. As with FOM, it uses the Arnoldi iteration to generate 
an orthogonal basis for the Krylov subspace. Since GMRES is a Krylov method, any 
vector X in the subspace can be written as 

X ^ xo + VmV (5.45) 

where xq is an approximate initial guess and x is the approximate solution to the 
system of linear equations. A residual vector is a measure of the accuracy of the 
approximate solution vector for a system of linear equations. The residual vector is 



50 



Table 5.6. FOM Algorithm 



1. Compute the residual vector Tq, with P := \ \ro\\2 and 
v^ 

2. Define the mx m Hessenberg matrix — hij=i^rn 
and initialize it to zero. 

3. For j = 1, m, Do 

4. Compute Wj := Mvj 

5. For i = Do 

6. hij = {wj, Vi) 

7. Wj = Wj - hijVi 

8. Enddo 

9. Compute /ij+ij = ||wj||2. If /ij+ij = 0, set m = j 
and compute the solution vector. 

10. Compute Vj+i — f^-^ 

11. Enddo 

12. Compute = if~^(/3ei) and = + VmVm- 



defined to be r = 6 — Mx where h is the right-hand side vector. The residual norm is 
the two-norm of the residual vector. It can be written as 

||r||2=|| \\2=\\h-M{xo + V^y) (5.46) 

Using the definition of the residual vector and Arnoldi iteration we can write 

r = ro - MVmV 

= (3vi-V^+iHmy (5.47) 
= Vm+i{f3ei- HmV)- 

Recall that in the discussion of orthogonal rotations, V is an orthonormal matrix. 
The residual norm is then 

||r||2 = \\b-M{xo + V^y\\2 (5.48) 

= \\Pe,-Hmy\\2. (5.49) 

The solution that GMRES produces, x, minimizes the residual norm. This can be 
found by finding the minimum residual solution with the vector ym- This vector is 
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Table 5.7. ALGORITHM :: GMRES(m) 



1. Compute tq = b — Mxo, P =|| tq ||2 and vi = tq/ j3 

2. For j = 1,2, ..,m, Do 

3. Compute Wj — Mvj 

4. For i = I, ...,j,Do 

5. hi J = {wj,Vi) 

6. Wj = Wj - hijVi 

7. Enddo 

8. hj^ij =11 Wj II2 If hj^ij = set m — j and goto step 11. 

9. Vj+i = Wj/hj+ij 

10. Enddo 

11. Define the (m + l)xm Hessenberg matrix Hm = {hij) i<i<m+i,i<j<m 

12. Compute ym, the minimizer of || /3ei — HmV II2, and = Xq + VmU', 



the minimizer of the residual norm in equation ( 5.46). The minimizer 

Um = min II /?ei - H^nV h (5.50) 

is computed by an inexpensive (m+1) x m least-squares problem, m is small for a 
practical LQCD application. The Arnoldi iteration used above minimizes the solution 
vector of the system of linear equations [56]. This gives the GMRES(m) algorithm 
listed in Talble 5.7. In the GMRES(m) algorithm, Givens rotations are employed in 
practice to determine the matrix elements hij. 

5.7.2 Restarted GMRES 

In practice, the GMRES algorithm becomes impractical due to growth of mem- 
ory and computational resources when the dimension of the Krylov subspace becomes 
large. As the dimension m increases, the computational cost increases at least as 
0{m^n) per cycle because of the orthogonalization of the elements of H. The mem- 
ory cost increases as 0{mn). [56] A solution to eliminate the high computational cost 
is to restart the Arnoldi iteration. The restarted GMRES(m) algorithm is listed in 
Table [56] The algorithm has the ability to exit when the desired residual norm is 



52 



Table 5.8. ALGORITHM :: Restarted GMRES(m) 

1. Compute To = b — Mxq, j3 =|| ro II2 .and, Vi = tq/ (3 

2. Generate the Arnoldi basis and the matrix Hj^ using the Arnoldi algorithm 
starting with vi 

3. Compute which minimizes || /3ei — Hm II2, and x„j = xq-\- VmUm 

4. If satisfied then Stop, else set xq — Xm and go to Step 1. 



reached for a given subspace size. If the residual norm is not satisfactory, the old 
Krylov subspace is replaced with a new subspace generated with the restarted ini- 
tial guess Xm- The restarted GMRES(m) method is the basis for many algorithms. 
One variation of this algorithm that we will consider is a Restarted Deflated GMRES 
method. 

5.7.3 Deflated GMRES 

For large, sparse matrices new GMRES techniques are required when the matrix 
eigenvalue spectrum contains small eigenvalues. For example, the Wilson matrix in 
LQCD contains small eigenvalues that give rise to exceptional conflgurations and need 
to be addressed to give sensible results. Techniques have been developed to solve 
problems of this nature for LQCD [63,64]. One method is GMRES with Deflated 
Restarting. This is referred to as GMRES-DR(m,k) where m is the dimension of the 
subspace and k is the number deflated eigenvalues for the spectrum. 

5.7.4 An Invariant Krylov Subspace 

For GMRES to remain effective, augmentation of the subspace by Rayleigh-Ritz 
vectors should return a Krylov subspace as well. In this subsection we verify that the 
Krylov subspace is still a Krylov subspace under the augmentation of approximate 
eigenvectors [63] . Since we are using restarted methods, each pass through the Arnoldi 
iteration (equation 5.32) between restarts is referred to as a "cycle". It was shown 
by Sorensen [65] that if the implicitly restarted Arnoldi method is restarted with 
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approximate eigenvalues (Ritz vector), the new initial vector is a combination of the 
desired Ritz vectors that generated these eigenvalues. So the subspace 

Span{yi, 1/2, Vk, Vm+u Mvm+i: M^Vm+u M"'-^v^+i) (5.51) 

is the implicitly restarted Arnoldi space in [66]. The vector Vm+i is the last Arnoldi 
vector from the previously-generated Arnoldi cycle. This vector is now the starting 
vector for the newly restarted Arnoldi cycle. It can be shown that equation ( 5.51) is 
equivalent to 

Span{y^, y2, Vk, My^, M^, M^^-^^y^), l<i<k, (5.52) 

where we have used the Arnoldi iteration from equation ( 5.32). Equation ( 5.51) is 
a Krylov subspace generated by a Ritz vector for each cycle. Similarly, in a restarted 
GMRES method, let Tq be the residual vector from the previous cycle. Then, the 
subspace is 

Span{ro, Mro, MVq, M^^-'^-Vq, m, m, Vk), (5.53) 

where yi are harmonic Ritz vectors. As shown in [66, 67], this subspace is equivalent 
to a subspace with the Harmonic Ritz vectors at the front of the subspace 

Span{yi, y2, yk, Myi, M%, M'^-^yi), (5.54) 

for 1 < i < /c. The span of these vectors is a Krylov subspace including the harmonic 
Ritz vectors, jji. By the preceding arguments we find that a Krylov subspace is still 
a Krylov subspace under augmentation of approximate eigenvectors. 

GMRES-DR(m,k) begins with a cycle of GMRES (m) which computes the so- 
lution vector and the matrix Hm- When the first cycle is finished, k-Harmonic Ritz 
vectors have been computed along with the matrix Vm using the Arnoldi recurrence. 
Vm is constructed by the vectors that span the subspace in equation ( 5.53). For the 
second cycle of GMRES-DR(m,k), as seen in equation ( 5.54), the first k-columns 
of the new matrix Vk consist of the orthonormalizied harmonic Ritz vectors. The 
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vector Vk+i must be generated by orthogonalizing the residual vector from the first 
cycle with respect to the columns of V^. Now that we have all the vectors needed 
to use the Arnoldi iteration we can form the rest of the Krylov subspace in ( 5.54). 
The GMRES-DR(m,k) algorithm [63] is summarized in table (5.9). It is important 
to realize that after the first cycle the Arnoldi iteration has changed slightly. The 
matrix Hm is upper Hessenberg except for a full leading A; + 1 by A; + 1 portion from 
the augmented eigenvectors. 

Computationally, it is reasonable to generate Schur vectors instead of eigenvec- 
tors. It is known that for any square matrix M, there exists a unitary matrix Q such 
that 

Q^MQ = R. (5.55) 

The Schur decomposition is then 

MQk = QkRk, (5.56) 

where the matrix R — RiR2...Rk is triangular and similar to M. The matrix Q is 
Q — Q2, Qk}- This is known as the Schur decomposition of M. Recall that if R 
is triangular and similar to M, then the diagonal elements of R are the eigenvalues 
of M. The columns of Qk are the Schur vectors of M, and they will be used as the 
approximate eigenvectors in this algorithm [63]. 

A simple example is useful to see how deflation is beneficial to solving a system 
of linear equations. In this example, we will augment the Krylov subspace with one 
approximate eigenvector. Let the source vector b = PiZi + + ••• + Pn^n- After 
the first cycle of standard GMRES the subspace is 

K = Span{zi, ro, Mro, MVq, M'^ro} (5.57) 

where ro is the new starting vector from the first cycle and zi is an exact eigenvector. 
The residual vector of this second cycle is generated by this Krylov subspace. The 
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Table 5.9. ALGORITHM :: GMRES-DR(m,k) 



1. Start: Choose m, the maximum size of the subspace 
and the desired number of 

approximate eigenvectors. Choose an initial guess, Xq, and 
compute ro = b — Mxq 

The new problem is M{x — xq) — tq. Let Vi — ro/ \\ ro \\ 

and /9 =11 ro || 

2. First cycle: Apply standard GMRES(m): use the Arnoldi 
iteration to generate Vm+i and Hm- Then solve the small min. res. 

problem min \\ c — H,nd \\ 

for d, where c = (3ei. Form the new solution vector 

Xm = xo + Vmd. Let f3 = hm+i,m, Xq = Xm, and ro = b- MXra- 

Compute the smallest k eigenpairs {Oi,gi) 

of + f3'H-^e^el. 

3. Orthonormaliztion of first k vectors: 
Orthonormalize the Harmonic Ritz vectors, gi 
and form an m x k matrix Pk- 

4. Orthonormaliztion of k + 1 vectors: Append a zero entry to each vector 
in the matrix to make them length m + 1. Then orthonormalize the 
short residual vector, c — Hmd against all the vectors in P^ to form Pk+i- 

5. Form portions of new H and new V 
using old H and old V: Let 

= Pl+iHmPk and y,-- = K„+iPfe+i. Then 
let Hk = Br"" and Vk+i = V^^f. 

6. Reorthogonaliztion of k + 1 vector: Orthonogalize v^+i against the 
earlier columns of the new Vk+i matrix. 

7. Arnoldi iteration: Apply the Arnoldi iteration from this point to 
form the remaining columns of V^+i and Ilm- Let P = /im+i.m- 

8. Form the approximate solution: Let c = V^^_iro and 
solve min || c — Hmd \\ for d. Let the new solution vector 
be x,n = Xq + Vmd. Computc the residual vector 

r = b- MXm = Vm+i{c - Hmd)- 

Check II r || = || c — H^d \\ for convergence, 

and proceed if not satisfied. 

9. Eigenvalue computations: Compute the k smallest eigenpair 
{ei,gi) of Hm + P'^H-'^emel. 

10. Restart: Let xo — Xm and ro = r. Proceed to Step 3. 
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solution vector that is spanned by this space is x — ^zi + g(M)ro where 7 is a free 
parameter. The residual vector for this cycle is 

r = {piZi + ...+l3nZn)-M{jZi + q{M)ro) 

= ((3iZi + ... + (3nZn)-^azi-Mq{M)ro (5.58) 
= {(3i- Xi^zi)- Mq{M)ro 

Notice that if we choose 7 = Pi/Xi the zi component does not contribute to the 
residual vector. By this choice of 7 the polynomial can "focus" on the rest of the 
spectrum. This approach is an alternative method to similar algorithms that use 
matrix preconditioning built of approximate eigenvectors to speed up the convergence 
of the residual vector [68-71]. 

5.7.5 Lanczos Method 

Krylov subspace methods rely on some form of orthogonalization of the Krylov 
vectors in order to compute an approximate solution to a system of equations. An- 
other class of Krylov methods are based on a biortliogonalization of a set of basis 
vectors. These projection methods are not orthogonal. The algorithm proposed by 
Lanczos [56] for non-symmetric matrices builds a pair of bases for the two subspaces 

Krn{M, vi) = span{vi, Mvi, M'^-^vi) (5.59) 

and 

K^{M^, wi) = span{wi, M^wi, {M^)"'-^Wi). (5.60) 

The pair of bases are built by the algorithm in Table 5.10. The scalars Sj+i 
and Pj+i are scaling factors for the bases vectors Wj+i and Vj+i respectively. If these 
scalar values tend toward a zero value in Steps 7 and 8, the algorithm will cease to 
converge and should exit in line 7 of the algorithm. 

As a result of lines 9 and 10, it is necessary to impose the constraint that 

= {vj+i,Wj+i). (5.61) 



57 



Table 5.10. The Lanczos Biorthogonalization Procedure 



1. 


Choose two vectors vi and wi that are parallel such that {vi,wi) = 1 


2. 


Set /?i = 5i = 1, Wo = Vo ~ 0. 


3. 


For j — 1, m, Do 


4. 


= iMvj,Wj) 


5. 


Vj+i = Mvj - ajVj - pjVj-i 


6. 


hatwj+i = M'^Wj — OijWj — SjWj-i 


7. 


Sj+i = \ {vj+i,Wj+i)\2. If Sj+i pa Stop. 


8. 




9. 




10. 





If equation (5.61) is satisfied we can write the tridiagonal matrix 



S2 0(2 Ps 



\ 



^m— 1 ^m— 1 /^r) 



V 



J 



Notice that the 5j is determined by the two norm of Vj and Wj and therefore are 
always positive. The (3j scaling parameter is then ±5j. 

It has been shown that if the Lanczos Biorthogonalization algorithm has not 
broken down by the m*^ step and {vi}^^^ m ^ basis of Km{M,vi) and {wi}^^^ ^ 
is a basis of K^{M'^ ,vi), then the following equations hold [56] 



W^MVm 



T 



(5.62) 
(5.63) 
(5.64) 



The and matrices can be interpreted as the projection matrices of M and M-^ 
onto the subspace K^{M, vi) and its orthogonal space K^{M'^ , vi). In practice, there 
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are many techniques which do not use the matrix M'^ , thus reducing the overhead of 
the Lanczos algorithm. These are referred to as transpose-free methods. 

5. 7. 6 Biconjugate Gradient Method 

The Biconjugate Gradient method (BiCG) is a non-symmetric Lanczos method. 
Imphcitly, BiCG solves not only the original system of equations, Mx — 6, but also the 
dual linear system of equations M'^x* = b*. The vectors wi and fi are not orthogonal 
to each other such that {vi,wi) ^ 0. wi is obtained from the initial residual vector 
b* — M^x*. The approximate solution vector that is obtained from the BiCG method 
has the form x = Vmdm, where dm = ^m^(/5ei). To find the inverse of the tridiagonal 
matrix we employ a LU factorization giving 

= {LmUj-\ (5.65) 

Now define the matrix = V^Uz} such that the solution vector is written 



Xm = Xo + VmT-\Pe{) (5.66) 
= Xo + PmLm'{Pe,). (5.67) 

The residual vectors for both the linear system of equations and its dual are denoted 
by rj and and are in the same direction as Vj+i and Wj+i respectively. 
For the dual system we define the matrix 

p; = ly^L-^. (5.68) 

Using the LU factorizations and the definitions in equations 5.68 and we can show 

{PmfMPm = L-^WZMVmU-' (5.69) 
= L-^T^C/-^ (5.70) 
= I. (5.71) 
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Table 5.11. ALGORITHM :: BiConjugate Gradient Method 



1. Compute the residual vector, Vg = b — Mxo and 
choose the dual residual such that (ro, r* 7^ 0) 

2. Set po = To and p* = r* 

3. Do J = 0, 1, ... convergence 

4. a, = (r„r*)/(Mp„p*) 

5. Xj+i = Xj + 

6. Tj+i = — ajMpj 

7. r*+i = r* - ajMp* 

8. /3, = (r,+i,r*+i)/(r,-,r*) 

9. pj+i = rj+i + (5jPj 

10. p*+i=r*+i + /?,p*. 

11. Enddo 



When equation 5.69 is true we say that the columns of and Pm are M-conjugate. 
We now have all the pieces to construct the BiCG algorithm for the system of equa- 
tions of M. This algorithm is found in Table 5.11. To solve the dual system, the 
residual vector in line 1 is replaced by r* — h* — M^x* and the update to the solution 
vector in hue 5 is x*j^-^ — x* -\- ajp*. 

5.8 GMRES Projection Method for Multiple Right Hand Sides 
In many physical applications, including lattice QCD, it is desirable to solve 
the same matrix equation for multiple right hand sides. It is important to solve these 
systems of equations together and take advantage of the fact that each right hand 
side shares the same matrix. We next describe the GMRES Projection method for 
Multiple Right-Hand Sides (RHS). 

This GMRES method employees GMRES-DR(m,k) to solve the initial system 
of equations (first RHS) and then uses the eigenvector information for the subsequent 
right-hand sides [72]. Specifically, projections over the approximate eigenvectors gen- 
erated in the GMRES-DR(m,k) algorithm are alternated with cycles of GMRES(m). 
This method is called GMRES(m)-Proj(k) where m is the dimension of the Krylov 
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Table 5.12. ALGORITHM :: MinRes Projection for Vk 



1. Let the current approximate solution vector be Xo with the associated 
system of equations M{x — Xo) = To- Let V^+i and 

come from equation 5.72. 

2. Solve the least squares problem min||c — Hkd\\i where 

3. Form the new approximate solution vector x = Xo + Vkd. 

4. Form the associated residual vector r — Vo — MVkd = — Vk+iHkd. 



subspace and k is the number of approximate eigenvectors used in the projection 
cycles. 

The approximate eigenvectors from GMRES-DR(m,k) span a small Krylov sub- 
space. These eigenvectors are generated by the "Arnoldi-like" recurrence 

MVk = Vk+iHk, (5.72) 

where is a n x k orthonormal matrix, Vk+i is similar to Vk with the inclusion of 
one extra row, and Hk is a full k+l x k matrix. The columns of Vk span a Krylov 
subspace as well as the subspace of approximate eigenvectors. 

A projection method is required over the approximate eigenvectors for the 
GMRES(m)-Proj(k) algorithm. A MinRes Projection projection over the subspace 
spanned by the columns of T4 is presented below (see Table 5.12). This projection 
is relatively inexpensive requiring 3/c + 2 vector operations (dot products and vector 
additions) of length n. 

The GMRES(m)-Proj(k) method is applied to all right hand sides except for the 
initial system of equations that is solved with GMRES-DR(m,k). The GMRES(m)- 
Proj(k) algorithm is presented in [72]. GMRES(m)-Proj(k) is summarized in Table 
5.13 The superscript on the solution and residual vectors indicates the current right 
hand side that is being solved. Notice that the projection in Step 3 adds very lit- 
tle overhead to the overall GMRES(m)-Proj(k) algorithm. A cycle of GMRES(m) 
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Table 5.13. ALGORITHM :: GMRES(m) - Proj(k) 

1. For the i*'* system of equations let M(x^^^ — x^^^°) = r*. 

2. If the right hand sides are related, project over the previous solution 
vectors. 

3. Apply the MinRes projection for V^. 

4. Apply one cycle of GMRES(m). 

5. Test the convergence of the current residual norm. If not satisfied, go to 
step 3. 



requires + 2m length n vector operations as well as the cost of m matrix- vector 
products. Step 3 requires 3k + 2 vector operations with no matrix- vector products. 

5.9 Multiple Shift Krylov Methods 
Recall our system of equations, 

Mx = b. (5.73) 

Using a Krylov subspace as a basis has many advantages for solving linear equations. 
One of these advantages is that it allows for multiple shifted systems of equations of 
the form 

(M - ai)x = r„ (5.74) 

to be solved simultaneously. There are more than three well-known methods that 
have been developed to solve this type of problem. This section is a review of Krylov 
methods for multiply-shifted problems. 

Let the Krylov subspace of dimension m generated by M and b be 

K^{M, b) = span {b, Mb, M"'-^b} . (5.75) 

The Krylov subspace is invariant under the shifts Uj, 

Km{M, b) = Krn{M - Ui, b)i = 1, US, (5.76) 

where ns is the number of shifted systems. 
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Krylov methods are iterative techniques for solving equation 5.73, where the m*'* 
iterate Xm satisfies Xm = Xg + Km{M, To) with = b — Mxq- If Xq = 0, then = b 
and we conclude from equation 5.76 that solution vectors of the shifted system of 
equations can be obtained from the same Krylov subspace as the unshifted system [73] . 
An important result of the shifted methods is that there is no added overhead to the 
calculation by adding the shifted systems, since the original subspace can be employed. 
Therefore, we can simultaneously solve the shifted system of equations for free. 

It was shown in reference [73] that to solve simultaneous shifted systems, resid- 
ual vectors of each system have to be coUinear such that r shift — ctrj, where is 
the initial system and a is a constant to be determined. Once a is determined, the 
residual of the shifted system is obviously a multiple of the initial residual. Therefore, 
the residuals come from the same Krylov subspace. Effectively, this corresponds to 
keeping the subspaces K{M,b) and K{M — ai,b) identical. 

5.9.1 Shifted BiCG 

BiCG employees a coupled two-term recurrence formalism. The two-term re- 
currence computes the polynomials p^, solution vectors Xm, and their residuals 
as seen in algorithm 5.11. If we demand that the polynomials of the shifted and 
unshifted system of equations are parallel, we must insist that the scalars (3 and a 
are coUinear such that (3shift — Ci(3 and aghift — C2OL. Ci and C2 are determined via 
the constraint that the two-term recurrence polynomials are coUinear for all shifted 
systems of equations. The details of how to determine the coUinearity coefficients for 
shifted BiCG are in [73]. The algorithm is summarized in Table 5.14 This Lanczos 
method has the advantage of relying upon short recurrence relationships. However, 
when the matrix is non-Hermitian, the computation of each basis vector used in the 
method requires a multiplication with M and which results in added overhead to 
the algorithm. 
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Table 5.14. ALGORITHM :: Shifted BiCG 



1. Initialize all solution vectors x = x shift = and r = b. 

2. Do A; = 1,2,... 

3. Employee the BiCG algorithm for the unshifted system of equations. 

4. For the shifted system, calculate Ci and C2 for Pshift and 
ashift, respectively. 

5. Determine the shifted polynomial Pshift = Pk^{<^)r - PshiftPshift 
where Pfc^(cr) is determined by the two term recurrence. 

6. Determine x^hift = x shift + ocsUftPshift 

7. Update the residual norms, r — r — aMp and r* — r* — aMp*. 

8. EndDo 



5.9.2 Shifted FOM 

The standard FOM method uses the Arnoldi algorithm to construct a Krylov 
subspace. Recall that and Hj^ are formed so that the first column of is 
Vi — ro/\\ro\\, and that the Arnoldi iterate is 

MVm = VmHm + hm+l,mVm+ie^- (5.77) 

For the shifted system of equations, the Arnoldi algorithm differs by a shift Ui as in 
equation 5.38 such that 

(M - (TiI)Vm = Vm{Hm " (^Im) + Wl,m^^m+ie^, (5.78) 

where 1^ is the identity matrix of dimension m. According to equation 5.78, the 
only modification to the original FOM method is that the small solution vector is 
solved via the system of equations {Hm — '^ImjUm — P^i [74]. 

It was also shown in reference [74] that shifted FOM can be restarted because 
Tjn is a multiple of the basis vector Vm+i- (ie = —hm+i,mVm+i{ym)m) The restarted 
shifted FOM algorithm is below. Similar to Shifted BiCG, restarted Shifted FOM only 
generates one basis {vi, ...,Vm} for all shifted systems to be solved simultaneously. 
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Table 5.15. ALGORITHM :: Restarted Shifted FOM 

1. Set To = b, (3l^= \ \ro\\, and a;^ = Xo- Set Vi = To/ Pl^- 
i is the current shift index. 

2. Construct Vm and for K^{M,vi). 

3. For all shifts construct 
Vm = (Hm - ailmY^eiPl^ 

Update the solution vector x'\^ = x\^ + VmVm- 

4. Exit if last shift has been computed to convergence criteria. 

5. Set (51^ = -hm+i,m{y\n)m for each shift. 

6. Set vi <— Vm+i to restart. 

7. EndDo 



5.9.3 Shifted GMRES 

Similar to both shifted BiCG and FOM, the same basis vectors for Km{M,b) 
and K^{M,b) can be used for the initial and shifted systems for GMRES so that 
the systems can be solved simultaneously where for convience we have defined M — 
(M — a I). However, for the restarted shifted GMRES method, rj and fj may not be 
parallel. This differs from the previous methods in that the matrix multiplications 
can not be saved upon the restart of the algorithm. A solution to this problem is 
presented in [75]. 

Any vector from the affine Krylov subspace can be written 

Xm^Xo+ Pm-i{M)ro, (5.79) 

where Pm-i is a polynomial of degree < m — 1. The corresponding residual is = 
b — Mxm- The residual can also be written in terms of the polynomial p^_i(M): 

rm = ro - Mp^_i{M)ro, (5.80) 
= gm(M)r«, (5.81) 

(5.82) 

where qm{M) = I — Mpm-i{M) with the initial condition that g(0) = 1. Similarly, 
we can define the residual norm and solution vector of the shifted system in terms of 
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this polynomial, 

Xm = Xo + qm-\{M)ro, (5.83) 
Tm = qm{M)fo, (5.84) 

where q is similar to the polynomial above but the identity has been shifted by cr. If 
we assume that the initial residual vectors are coUinear, then fo — oloTq where e C. 
With this initial condition, the constraints to keep the shifted systems parallel are 

Tm = ttm'^m, (5.85) 

which yields 

OLoqm{M)ro = amqmiM - aI)ro. (5.86) 

Equations (5.85) and (5.86) are the defining equations for q^ and am- 

The Arnoldi equation for the initial and shifted system of equations are 

MVm = Vm+lBm (5.87) 

and 

MVm = Vm+lHm, (5.88) 

with Hjn = Hm — ailm+i,m where the last row of the m+1 x m identity matrix is full 
of zeros. 

The coUinearity condition for the residual norms in equation 5.85, along with the 
definition of the shifted Arnoldi recurrence in equation 5.88 yield the underdetermined 
equation [75] 

HrnVm + (/3ei - Hmym)oim = Qiol |ro| ^Bl. (5.89) 

There are two unknown variables, ym and a^, in equation 5.89. In practice, a QR 
factorization can be used to solve this equation to determine a^.- Once the coUinearity 

parameter is determined, the solution vector {jm is determined, which allows the 
shifted GMRES algorithm to be restarted. 
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Table 5.16. ALGORITHM :: Shifted GMRES(m) 



1. Set the initial guesses for Xo, Xo such that fo = Oo^o- 

2. Employee the Arnoldi algorithm to determine Vm and Hj^. 

3. Use GMRES(m) for the initial system of equations to determine 
ym, the sohition Xo + VmUm and 

m I ■ 

4. Using the output of the initial GMRES(m) algorithm, 
solve equation 5.89 for and ym- 

5. Determine the solution of the shifted system, x^ — Xo + Vmym- 



Table 5.17. ALGORITHM :: Restarted Shifted GMRES(m) 



1. Set the initial guesses for 

the shifted GMRES(m) algorithm. Set the restart value, k. 

2. Do j = 1,2, 

3. Use Shifted GMRES(m) to determine 
and oP-^^ such that f]^^ = ai^^rl^'^ 
via equation 5.89. 

4. Stop if the residual norm of the initial and shifted systems 
has reached convergence criteria. 



For a fixed m value in GMRES(m), the Shifted GMRES(m) algorithm is writ- 
ten in Table 5.16 The restarted shifted GMRES(m) algorithm is in Table 5.17. An 
example comparing restarted, shifted GMRES(m) and other shifted methods is shown 
in the next section in Figures (5.2) and (5.3). 

5.10 Krylov Methods with Multiple Shifts and Multiple Right Hand Sides 
We have already discovered that GMRES-DR(m,k) can be implemented in ap- 
plications where it is desired to solve multiple shifts simultaneously or multiple right- 
hand sides. Sometimes both of these methods are needed for the same application. 
One specific application is Wilson LQGD. In the Wilson formahsm the different shifts 
represent different quark masses and the right hand sides are quark sources on the 
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lattice. A new deflated GMRES technique for both multiple RHS and multiple shifts 
is provided in this thesis. 

Consider the large system of linear equations that not only has multiple right 
hand sides, but also has multiple shifts for each right hand side. Let nrhs be the 
number of right hand sides and ns be the number of shifts. The problem to solve is 
then 

(M - ai)xi = bj, (5.90) 

with j — 1, nrhs and i — 1, ...,ns. M is a large matrix which may be nonsymmetric 
or complex non-Hermitian. In this section ai is referred to as the base shift. 

5.10.1 Shifted GMRES-DR(m,k) 

To achieve our goal of an algorithm to solve multiple shifts simultaneously for 
multiple RHS, another new shifted GMRES (m) algorithm was needed which would 
take advantage of the deflation of eigenvalue spectrum to speed convergence. Such 
a method has been developed in this work and is referred to as GMRES-DRS(m,k), 
which stands for shifted GMRES-DR(m,k) where m and k have the same meaning as 
in the GMRES-DR(m,k) algorithm. 

For the deflated GMRES (m) method wc can solve multiple shifted systems 
concurrently. The derivation is the same except that we employee the "Arnoldi- 
like" recurrence that differs from the Arnoldi recurrence in that Hjn is an upper 
Hessenberg m+1 x m matrix except for a full k + lxk leading portion that contains 
approximate eigenvector information needed by the deflation technique. The shifted 
GMRES-DR(m,k) algorithm is listed in Table 5.18. GMRES-DRS(m,k) is related to 
the method GMRES-E [76]. In practice, the QR factorization in this algorithm uses 
Givens rotations to solve the system of equations in Step 4. 

A comparative study between shifted GMRES(25) and the new GMRES-DRS(25,10) 
algorithm was carried out to study the convergence of each method. The matrix had 
n — 1000 and is bidiagonal with 0.1, 1, 2, 3, 998, 998 on the main diagonal and I's 
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Table 5.18. ALGORITHM :: Shifted GMRES-DRS(m,k) 



1. At the begging of a cycle of GMRES-DRS(m,k), assume the current 
problem is (M - ail){xi - Xo,i) = f^ir^i, 

with Pi — 1 and where Xg^i is the current 
approximate solution to the i*'' shifted system. 

2. Apply GMRES-DR(m,k) to M and generate the "Arnoldi-like" recurrence 

3. For the base system, solve the minimum residual problem 
min| |c — {Hm — ail), where 

c = V^+i'^o,! and /isam + lxm identity matrix. 

The new approximate solution vector is Xi = Xo,i + Vmdi- The new 

residual vector is ri = r^.i — MVmdi = Vo^i — Vm+iHmdi. 

4. For the other shifted systems i = 1, ...,ns form s = c — {Hm — <JiI)di. Apply 
a QR factorization: — ail — QR. Solve Rdi — PiQ^c + aiQ'^s, 

using the last row to solve for ctj and the first m rows for dj. 

5. The new approximate solution vector of the i*'* system is x = Xo,iVmdi, and 
the new residual is = aiVo- 

6. Test the residual norm for convergence. If not satisfied, for i = 2, ns set 
Pi — ai and for i — ns, set Xo^i — Xi and r^j = rj. 

Then go back to step 1. 



on the superdiagonal. Thus, M has the form: 



M 



.1 1 .. 

Oil:: 
2 1 

.. .. 



The associated right hand side is randomly generated. The shifts are a = 
1,-0.4,-2. This matrix has a small eigenvalue that slows down the convergence of 
the residual norm for restarted, shifted GMRES(m), especially for the base system. 
Restarted, shifted GMRES(m) is compared to GMRES-DRS(m,k). The results are 
in Figure (5.2). 
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Figure 5.2. Comparative results for restarted, shifted GMRES(m) and GMRES-DRS(m,k) 

Shifting the matrix by 0.4 improves the convergence of GMRES(25) because the 
smallest eigenvalue in the spectrum is moved from 0.1 to 0.5. GMRES-DRS(25,10) 
converges rapidly once it generates approximate eigenvectors corresponding to these 
eigenvalues. The convergence for all three shifted systems is similar because once the 
small eigenvalues are removed by the deflation, shifting has little effect on convergence. 

In this example the second and third shifted systems converge faster than the 
first base system. However, in some situations, there can be convergence problems 
for the non-base systems [74]. This author compare multiply-shifted GMRES(m) and 
FOM. Since FOM has parallel residuals for all shifted systems, it is argued that it is 
a more natural approach to the shifted problem [74]. However, convergence depends 
on the roots of the polynomial, which correspond to eigenvalues of M and where the 
roots fall in relation to the shift. The next example demonstrates this dependence. 
Recall that shifted GMRES(m) uses the same polynomial for all shifted systems. 
Likewise, shifted FOM also uses one polynomial for all shifts. This polynomial is 
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subject to the condition g(0) = 1. We may view this as the polynomial being 1 at 
zero and the spectrum shifted, or the polynomial has the value 1 at the shift and the 
spectrum fixed as that of M. We prefer the second view. So, for shifted GMRES(m) 
with the base system (M — o"j/), we view the polynomial chosen by this method as 
being 1 at ai and needing to be small over the spectrum of M. 

In the next two examples, plots are given of the roots of the polynomials. 
Because the polynomial must be small over the spectrum of M, a small root of the 
polynomial at the shift can cause a convergence problem. So for Krylov methods to 
be effective, the roots of the polynomial need to generally stay away from the shift 
value. 

For the bidiagonal matrix used in the first example, we apply shifted GM- 
RES(40) and shifted FOM(40) with shifts a = 0.4,0. The results are in Figure (5.3). 

For the base shift of a = 0.4, GMRES(40) works better than FOM(40). Alterna- 
tively, FOM performs better than GMRES for the shift value a = 0. Figure (5.4) 
shows the harmonic Ritz values nearest the shift for 25 cycles of GMRES. These are 
the roots of the GMRES polynomial for the system (M — 0.4/)a;i = b, shifted so they 
correspond to the spectrum of M. 




Figure 5.3. Comparative results for shifted GMRES(40) and shifted FOM(40) 
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Figure 5.4. Distribution of Smallest Harmonic and Regular Ritz values, m=40, 25 cycles. 

The harmonic Ritz values avoid the region about 0.4. This is a good result since 
the shifted GMRES(40) polynomial for a = 0.4 is defined to have a value of 1 near 
this value and be small over the spectrum of M. This polynomial cannot be effective if 
there is a root near 0.4. GMRES(40) converges slowly for this fairly difficult problem. 
In contrast to GMRES(40), FOM(40) with a = 0.4 is not able to converge because 
the roots of the FOM polynomial are not separated from 0.4 as seen in the bottom 
plot in figure (5.4). For the second shift value of cr = shifted GMRES(40) is not 
effective. There are too many harmonic Ritz values at and surrounding 0. Shifted 
FOM(40) gives erratic convergence because of some Ritz values near zero, but it does 
manage to converge. See [74] for more comparative examples between shifted FOM 
and GMRES. We next give an example where GMRES is more effective than FOM. 

The matrix is the same as in the previous two examples. The base shift is cr = 
and a second shift of cr = 1.4 is used. Shifted GMRES (80) is compared with shifted 
FOM(80). The results are shown in Figure (5.5). 

Some of the Ritz values in the FOM(80) spectrum fall around 1.4, while there is 
a gap around these values in the harmonic Ritz values. These approximate eigenvalue 
distributions are shown in Figure (5.6). 
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Figure 5.5. Comparative results for shifted GMRES(80) and shifted FOM(80) 



Figure 5.6. Distribution of Smahest Harmonic and Regular Ritz values, m=80, 25 cycles. 

In Figure (5.6) we see that shifted GMRES(80) works much better for the 
second system. Therefore, we conclude that methods for shifted GMRES(m) and 
shifted FOM must both be used with caution. However, deflating eigenvalues can 
help this problem. Figure (5.5) also has a plot of FOM-DRS(80,2) for a = 1.4 and we 
see that deflating only two approximate eigenvalues fixes the convergence problem. 
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Deflated methods in general will eliminate the "small eigenvalue" problem. For the 
rest of the examples in this thesis the matrices and shifts are such that the base 
systems have the slowest convergence. 

5.11 Deflated GMRES for Multiple Right-Hand Sides and Multiple Shifts 
We now consider solving multiply-shifted systems that also have multiple right- 
hand sides. It is important to reuse information or share information among the 
right-hand sides. Block methods [56] share information between right-hand sides. It 
is possible to design multi-shifted versions of both Block-GMRES [56] and Block- 
GMRES-DR [77]. However, here we will concentrate on a non-block approach. The 
right-hand sides are solved separately, and eigenvector information from the solution 
of the flrst right-hand side is used to assist the subsequent ones. More speciflcally, we 
will generalize for multiple shifts the GMRES-Proj approach mentioned in Section 
5.8. See [72] for more on this method, including comparison with block methods. 

First, some of the difficulties of deflating for subsequent right-hand sides will be 
discussed. Suppose the first right-hand side has been solved and approximate eigen- 
vectors have been generated. Then for the non-shifted case, there are several ways 
to deflate eigenvalues. Some of these are given in [72]. However, generally they do 
not work for multiply-shifted systems. For example, if the deflation involves building 
a preconditioner from the approximate eigenvectors [72,78,79], then as mentioned 
earher, the differently shifted preconditioned systems cannot be solved together. 

For the GMRES (m)-Proj(k) method, there is trouble with one of the two phases. 
We know the GMRES portion can be adapted to keep right-hand sides parallel for 
multiple shifts. However, the phase with projection over approximate eigenvectors 
generally fails to produce parallel residual vectors. Even though this projection is 
over a Krylov subspace of dimension k spanned by the columns of V^, this subspace 
does not contain the current right-hand side (the residual vector); so the derivation 
in 5.8 does not work. Specifically, the residual vectors between adjacent right-hand 
side vectors are not parallel since ro,i ior i — 1, ...nrhs and ro,i are not in the span of 
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the columns of V^+i. One case in which the projection does keep the residual vectors 
parallel is exact eigenvectors, as shown below. 

Theorem :: Assume that before the minres projection, the shifted systems are 

M{x - Xi) = ro,i 

with ro,j = Aro,! for i = 2, . . . , ns. Then after the minres projection over the subspace 
Span{zi, Z2, ■ ■ ■ , Zk}, where zi through Zk are eigenvectors of M, the residual vectors 
are still parallel. Proof. Let Zk be the matrix with zi,...,Zk as columns. For 
exact eigenvectors, the minres projection is equivalent to Galerkin. The Galerkin 
orthogonality condition gives for the i*'*-system 

Z^{(3iro,i-{M-aiI)Zkdi) = 0. 

Solving gives 

di = A(Afc - aih)-\ZlZk)-'Zlro,i, 

where is the k x k diagonal matrix with diagonal entries Ai through and is 
the k X k identity matrix. The residual vector after projecting is then 

<r = Piro,i-{M-aiI)Zkd, 

= Aro,i - Pi{M - aJ)Zk{Kk - aJkY^Zl ZkY^ Zlr^,^ 

= l3iro,i - l3iZk{Kk - aiIk){Ak - aiIky^{Zl Zk)'^ ZItq^i 

= A(7 - Zk{ZlZk)-^Zl)ro,,. 

This shows that all ro,i are multiples of each other. 

So, one option for using GMRES-Proj with multiple shifts is to use only fairly 

accurate eigenvectors. We could sort through the approximate eigenvectors computed 
by GMRES-DR and apply only ones with acceptable accuracy to the projection in 
GMRES-Proj. This is now tested. 

Example 1. For the same matrix, wc first solve the ai = system to an accuracy 
relative to the residual norm below 10~^° or 250 matrix-vector products. Table 5.19 
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Table 5.19. Effect of projecting over different accuracies of eigenvectors 
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8 
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3.7e-5 


7.1e-6 
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1.3e-9 


6 


8.0e-5 
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2.8e-6 


3.3e-10 


180 


1.3e-9 


4 


l.Oc-6 


255 


5.4e-8 


2.4C-11 


255 


5.5C-10 


2 


5.4e-9 


435 


1.5e-8 


9.9e-12 


435 


3.6e-10 



shows the accuracy of the eigenvectors thus produced in its second column. The worst 
residual norm of the 10 approximate eigenvectors is 3.3 x 10~^. If, for example, only 
the four smallest approximate eigenvectors are used, their accuracy has a residual 
norm of 1.0 x 10~^ or better. Table 5.19 lists the number of matrix-vector products 
for the relative residual norm of the base shifted system with the second right-hand 
side to reach 1.0 x 10~^. We see that convergence is better using all ten eigenvectors. 
However, the fourth column gives the accuracy attained by the worst of the last 
two shifted systems (usually it is the third shift). It only reaches a residual norm 
of 3.3 X 10~^ if all 10 approximate eigenvectors are used in the projection. With 
only four eigenvectors, the residual norm reaches a better level of 5.4 x 10""^, but the 
convergence is almost twice as slow. 

The problem with this approach is that the eigenvector computation during 
solution of the first right-hand side needs to be done to a considerable accuracy, since 
we do not want to slow down convergence of the subsequent systems. If many right- 
hand sides are to be solved, this extra expense might not be significant. However, we 
next propose an approach that does not require eigenvectors to be accurate. 

The key idea is that although the residual vectors cannot be kept parallel, they 
can be chosen so that they relate to each other. We force the residuals of the non-base 
systems to be parallel to the residual of the base system except for a component in 
the direction of Vk+i, the last column of the 14+1 matrix from Equation (5.72). 

Now we will look at the correction phase that is needed at the end of GMRES- 
Proj. Assume that we have already solved shifted systems with the extra right-hand 
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side Vk+i (this solution will be discussed next) and have 

(M - ail)si = Vk+i. (5.91) 

We assume that for a particular right hand side, the systems have been solved by 
GMRES-Proj to the point that the residual is only in the direction of Vk+i'- 

h-{M - ail)xi = ri = 'jiVk+i, (5.92) 

for i = 2, . . . , ns and for some scalar 7j. Here Xi is the approximate solution to Xi 
(the superscripts for the particular right-hand sides are left off here for simplicity). 
We perform a Galerkin projection for system (5.92) over the subspace spanned by the 
single vector Sj from solution of (5.91): 

sJ{M - (Til)si5 = -fisjvk+i- 
Using Equation (5.91), this becomes 

Then 5 — ^i. To determine 7j, we start with r = ^iVk+i- Multiplying both sides by 
v^j^-^ and using that v^+i is of unit length gives 

7i = vl^ir. (5.93) 

The corrected solution of the system is Xi — Xi -\- 5si. 

We need to fill in the method for solving (5.91), the shifted systems with the 
extra right-hand side. First, GMRES-Proj is applied until the residual is negligible 
except in the direction of Vk+i- So we assume that 

(M - ail)si = r (5.94) 

is the current system, where 

r = Vk+i - (M - ail)si = -fiVk+i, 
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for some scalar 7j. Rearranging gives 

(M - ail)si = (1 - -fi)vk+i. (5.95) 

Applying Galerin projection over the subspace spanned by the single vector Si to the 
system (5.94) gives 

sJ{M - ail)si5 = -fisjvk+i. 
With Equation (5.95), this becomes 

(1 - 'yi)sJvk+iS = 'fisjvk+u 

and this simplifies to 

1 -7i 

So the corrected solution is Sj = Sj + Y3^Si- Now 

1 

Si ~ Si. 

1 -7i 

Finally, the ji is determined to be = vj+ir as it was for (5.93). 

We next give the algorithms for solution of the systems with second and subse- 
quent right-hand sides and for the extra right-hand side. Note these are in order of 
how they were derived here, not in order of how they are actually used. The algorithm 
for solution of the systems with second and subsequent right-hand sides is given in 
Table 5.20. 

The algorithm for solution of the systems with extra right-hand sides is given 
in Table 5.21. 

Example 2. We use the same test matrix. All right-hand sides are generated 
randomly. The systems with the first right-hand sides are solved with GMRES- 
DRS(25,10) as before. Then the extra right-hand Vk+i systems are solved (for all 
shifts) with GMRES(15)-Pr(10)-Sh. Finally, the second right-hand side systems 
are also solved with GMRES(15)-Pr(10)-Sh. All relative residual tolerances are 
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Table 5.20. ALGORITHM :: GMRES-Proj-Sh for the second 
and subsequent right-hand sides 



1. Consider the systems with the j'th right-hand side (and with all ns shifts). 
At the beginning of a cycle of GMRES(m)-Proj(k)-Sh, assume the current 
problem is (M - ail){xl - J = Aro,i, 

with Pi — 1, and where xo,i is the current 
approximate solution to the ith shifted system 

2. Apply the Minres Projection for to the first right-hand side. 
This uses the V^+i and Hk matrices developed while solving 
the first right-hand side with GMRES-DRS(m,k). 

3. For shifted systems is — 2 . . . ns, solve (iffc — ail)di = /3i{Hk — ail)di. 

4. Apply one cycle of GMRES(m)-Sh. 

5. Test the residual norms for convergence (can also test during GMRES cycles). 
For the non-base systems, can ignore the error term in the direction of v^+i- 

6. Correction phase: Suppose the computed solution for the ith shifted system 
Let the solution to the system with the extra right hand side 

so far is xj. Vk+i and shift (Tj be Sj. 

The corrected solution is (^xf)'^'^^''*^''' — xl + (t'^+i * r) * xvi. 
The corrected residual norm can now be calculated. 



Table 5.21. ALGORITHM :: GMRES-Proj-Sh for the extra right-hand side Vk+i 



Same as for previous algorithm except for . . . 

1. Consider the systems with right-hand side Vk+i (and with all ns shifts). 

2. Correction phase: Suppose the computed solution for the ith shifted system 
so far is Sj. The corrected solution is 

Si = Si + (j^) * Si, with 7j = v^^^r. 



rtol — 1.0 X 10~^. Figure (5.7) has residual curves for only two shifts, the base 
shifts of zero and a — —2. The dotted fine shows the uncorrected residual norm 
for the second shift, while the dash-dot line has the second shift residuals if they 
are corrected (actually the correction needs to be done only once at the end of each 
right-hand side). The uncorrected residual norm for the second shifted system lev- 
els off at 4.0 X 10~^, but this is fixed by the correction phase. The convergence is 
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faster than for GMRES-DRS, because the eigenvectors are used from the beginning to 
speed up the convergence. Also the cost of GMRES(15)-Proj(10)-Sh is less than for 
GMRES-DRS(25,10), because it is fairly inexpensive to project over the approximate 
eigenvectors compared to keeping the eigenvectors in the GMRES-DR subspace. Here 
the expense for the extra right-hand side is fairly significant, however it will not be 
if there are more right-hand sides. Figure (5.8) has the case of solving a total of 10 
right-hand sides. Also, the extra right-hand side is solved only to relative residual 
tolerance of 1.0 x 10^^. Now the expense for the extra right-hand side Vk+i is small 
compared to the amount saved by speeding up the solution for all the remaining 
right-hand sides. 




Figure 5.7. Solution of first rhs, extra rhs and second rhs with two shifts. 

Example 3. At the end of the previous example, the extra right-hand side is 
solved to low accuracy, but the correction for the subsequent right-hand sides is still 
successful. We now experiment with solving the extra right-hand side to different 
levels of accuracy. Table 5.22 shows the accuracy after correction for the a = —2 
system when the extra right-hand side system is solved to relative residual tolerances 
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Figure 5.8. Solution of first rhs, extra rhs and ten rhs witli two sliifts. 

ranging from 1.0 x 10~^ down to 1.0 x 10^^ (the tolerance is checked for termination 
only at the end of GMRES cycles). The first and second right-hand side systems 
are solved to three different residual norm tolerances (1.0 x 10-^ 1.0 x lO"® and 
1.0 X 10^^°) in the three rows of the table. The conclusion of this experiment is 
that the extra right-hand side systems do not need to be solved very accurately. 
With tolerances 1.0 x 10~^ and 1.0 x 10~® for the first and second right-hand sides, 
the extra right-hand side systems need only to be solved to a relative tolerance of 
1.0 X 10"'^ for essentially full accuracy. 

Example 4- We look at a Wilson-Dirac matrix from lattice QCD. The dimension 
is 393,216 by 393,216. The value of k is 0.158 for the base shift. This is approximately 
i^criticai- The right-hand sides are unit vectors associated with particular space-time, 
Dirac and color coordinates. Often there are a dozen or more right-hand sides asso- 
ciated with each matrix and perhaps seven shifts for each right-hand side. We will 
just show solutions of the second right-hand side for three shifts, a = 0, —0.3, —0.5. 
The first right-hand side is solved with GMRES-DRS (50,30) to a residual tolerance 
of l.e-8 and the extra right-hand side to 1.0 x 10^^. Then for the second right-hand 
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Table 5.22. Effect of solving the extra right-hand side system to different accuracies 
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side, GMRES-Proj uses 30 approximate eigenvectors for the projection in between 
cycles of GMRES(20). See Figure (5.9) for the results. GMRES(20)-Proj(30)-Sh can 
converge in about one-tenth of the iterations needed for GMRES(20). To reach a 
residual norm of less than 10~^ for the toughest system with shift of zero takes 2680 
matrix- vector products for GMRES(20)-Sh and 280 for GMRES(20)-Proj(30)-Sh. 

5.12 Projection Methods for tmQCD 
To solve shifted system of equations simultaneously each shifted system of equa- 
tions must use the same Krylov subspace. In tmQCD at maximal twist, the matrix 
M changes for each shift due to necessary even-odd preconditioning of the problem. 
Since the matrix changes for every shifted system, the Krylov subspace used to solve 
the base system will not work for the following shifted systems. So, simultaneous 
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Figure 5.9. Solution of second RHS for large QCD matrix with three shifts. 



82 

shifted Krylov methods do not work for the tmQCD formahsm and the systems must 
be solved serially. Preliminary results of a new method that help convergence of the 
subsequent shifted systems using a projection over solutions are presented in this 
section. 

The solution vectors of the shifted system of equations (M — ail) for i = 
1, j — 1 can help the convergence of a shifted system (M — ajl) where j > i. To 
help the convergence of the j*^ system of equations a MinRes projection over the 
previous j — 1 solution vectors is used to create an initial guess x^j for the current 
system. In this method, we solve the most difficult system of equations last. In doing 
so, we take advantage of the projection over all the previous solution vectors. 

Again, the problem referenced in Example 2 is explored. The six shifted systems 
in Figure (5.10) correspond to the shifts cr^ = {-0.05, -0.04, -0.03, -0.02, -0.01,0} 
for i = 1, 6. 




Figure 5.10. Residual vector for serial shifted systems as a function of MVP's for GM- 
RES(30) 

The initial base system is solved with GMRES(30). It is obvious from Figure 
(5.10) that the convergence is helped for each successive RHS. The first system of 
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equations required approximately 120 matrix-vector products to form a residual of 
10^''. In contrast, the last system of equations (which is the most difficult shifted 
system) needed approximately 50 matrix- vector products. Next, the same example 
using the same shifts is repeated for GMRES-DR(30,10). 



Mulitple Siiifts 
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Figure 5.11. Residual vector for serial shifted systems as a function of MVP's for GMRES- 
DR(30,10) 

In Figure (5.11) the total number of matrix-vector products to solve all of the 
shifted systems is approximately 380. The total number of matrix- vector products 
(MVP) in Figure (5.10) is approximately 450. The saved MVP are a result of deflation 
and the projection over the previous solution vectors. In future work, this method 
will be extended for multiple RHS vectors using GMRES(m)-Proj(k). 



CHAPTER SIX 
Disconnected Sea Quarks 

Disconnected loop calculations have historically been a challenging problem for 
hadronic physics. Exact calculations of light quark matrix elements at each lattice 
point is extremely expensive computationally and currently not realistic with our 
current computer resources. An alternative to the exact calculation is to utilize an 
unbiased, stochastic estimate of the loops [6,80,81]. 

This technique uses noise theory to project out the loop operator expectation 
values. A continuing challenge with the noise methods is to reduce the variance of 
the calculation such that a stronger signal is acquired. By reducing the variance, 
the computational costs also decrease for the operator calculation. Higher order 
subtraction results are presented in this chapter as well as preliminary subtraction 
results for a twisted perturbative subtraction technique. 

6.1 Noise Theory 
The disconnected loops can be described by the systems of equations 

Mx = T] (6.1) 

where M is the quark matrix, x is the solution vector and r; is a noise vector that is 
used to project the disconnected signal. The constraints on the system of equations 
are 

< rii >= 0, < r]ir]j >= Sij, (6.2) 

where the average is over all noises used. Using these identities, a particular inverse 
matrix element M^^, can be determined by 

<r]jXi> = J]m-^ < 77,-77fe > (6.3) 

k 

= Mr^. (6.4) 
84 
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At this time it is instructive to review the basic principals of matrix inversion 
using noise theory [6,81]. The expressions for the expectation value and variance of 
a matrix in terms of a general noise vector are found in Ref. [81]. 

Let the average of the projection of one general noise vector onto the other be 

1 ^ 

Xmn= ^^Vmk-nlk^ (6-5) 
k=l 

for (m, n = 1, ...,£)), where D is the dimension of the matrix and {k — 1,...,N), 
where N is the number of noises used. The matrix Xmn is hermitian with expectation 
value < Xmn >= Smn as abovc. Using this notation, the variance of the measurement 
is defined to be 



V[TrQX] = <\J2QmnXmn-TrQ\''> (6.6) 

m,n 

= XI < l^nn- ir >< \Qnn\ >^ 

n 

E/ \ l2i |2 ^ / \2\ 

(< \Xmn\ iQmnl + QmnQnm < [Xmn) >)■ 

6.1.1 Real Z(2) Noise 

The Z{2) noise constraints are 

< >= < {Xmnf >= (6.7) 

for m ^ n. 

Allow the matrix, QX, to be real. Now consider equation (6.6), when we apply 
the constraints for Z{2) noise and notice < \Xnn — Ip >= 0. We may write the 
variance as 

VlTvQXreal] = ^ J](kmnr + ^mnqlJ- (6.8) 

Therefore, by equation (6.8), Z{2) noise has the lowest variance of any real 
noise. This implies that the variance for Z{2) noise is a result of the off-diagonal 
matrix elements. 
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6.1.2 General Z{N) Noise 

For general Z{N) noise (A^ > 3) we have a different set of constraints, 

for m ^ n. Similarly to the real case, we have that < |X„„ — Ip >= 0. Thus, the 
expression for the variance becomes 

V[TrQXz^^)] = ^ 5^ \qran?- (6.10) 

For a general matrix Q, the variance relationship of Z(2) and Z(N) is not fixed. 
The difference in the variances is due to that fact that the square of an equally 
weighted distribution, as is the case for Z (2), is not itself always uniformly distributed. 
In contrast, the square of the uniformly weighted Z[N) for > 3 is uniformly 
distributed. Even so, if the phases of q^n and are not correlated, the variances 
for Z{2) and Z{N) {N > 3) are approximately the same. For the operators that we 
calculate this appears to be the case [6]. 

6.2 Perturbative Subtraction 
Perturbative noise subtraction gives a computationally efficient and effective 
way to reduce the variance of disconnected operators by using noise theory meth- 
ods [6]. A review of the methodology in reference [6] is useful for our twisted mass 
formalism. 

The trace of a matrix is obviously invariant under addition of another traceless 
matrix. Given two matrices, Q and Q, the expectation values are related by 

<Tr{{Q-Q)X}>^<Tr{QX}>, (6.11) 

where Q is traceless. The variance, however, is not invariant under the addition of 
the traceless matrix Q: 

V[Tr{{Q - Q)X] =< I Y,i<lrnn - qmn)Xmn - TrQ\^ > . (6.12) 
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The variance is completely determined by the off-diagonal elements of Q and Q. The 
variance can be minimized if the off-diagonal elements of Q and Q are similar. It 
is important to find a Q matrix that is traceless because we only want to reduce 
the variance of Q not the diagonal elements that contribute to the disconnected loop 
expectation values. 

The Wilson matrix can be written as 

(*^-')" - JJJ^y ("3) 

where the capital indices are over space, color and Dirac indices (/, J = {x,a,a}). 
For Wilson fermions, the matrix elements Pjj are 

PiJ = E[(l - ^^)UA^)Scc,y-a, + (1 + 7,)Ul{x - a,5,,y+a,)]. (6.14) 

Expanding equation (6.13) in a geometric series in the hopping parameter k, one can 
write the perturbative Wilson quark matrix as 

M^ertiP) ^I + '^P + H^P^ + l^'P^ + ■- (6.15) 

To reduce the variance of the weak matrix elements expectation values, < rjjM^^rjk >, 
it is natural to choose Q to be the perturbative quark matrix MpertiP)- According to 
equation (6.12), the variance is calculated from the off-diagonal elements of Q. With 
the perturbative Wilson matrix as our choice of Q, we calculate the difference of the 
expectation values of the Wilson quark matrix, < rjjM^^rjk >, and the perturbative 
matrix < r]jMpglf{P)ikr]k > to reduce the variance of the disconnected loop operator. 
However, the perturbative quark matrix that we used to reduce the variance is, in 
fact, not traceless and therefore the diagonal elements of M^^ are changed by the 
subtraction of M~J^^. 

In general, the expectation value of an operator O is 

< ijjOijj >= -Tr{OM-^). (6.16) 
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The subtraction that changed the variance of the loop operators also changed the di- 
agonal elements of the M matrix. These values contribute to the vacuum expectation 
value and need to be added back in to get the full, unbiased answer. One way to 
calculate the diagonal elements is to explicitly construct all the gauge invariant paths 
that contribute to a given operator. 

Only closed loop, gauge invariant objects contribute to the trace in equation 
(6.16). The local operators require a perturbative correction starting at 4*^^ order and 
non-local, vector operators need corrections starting at S^'^ order in /t. (The vector 
operator also requires a correction at zeroth order in n. ) Another way to view the 
perturbative trace is that only closed path objects with an area A contribute to the 
trace in equation (6.16). A general picture of the vector and scalar operators is in 
Figure (6.1). The quark lines in this figure represent all possible quark propagators 
that add to the operators. 




Figure 6.1. General diagram of the quark line contributions to the flavor singlet scalar and 
vector operators. 



89 

The local scalar operators that contribute to the operator signal begin and end 
on the same lattice site. Vector operators are next-neighbor interactions that are 
connected with a gauge link U{x). 

6.2.1 Subtraction Methods 

To reduce the noise variance of the disconnected operators we need to subtract 
the inverse of the pcrturbativc quark matrix M^^^ from the quark matrix Mjj . In 
practice, we construct the perturbative matrix to a desired order in k with equation 
(6.15). 

The construction of the perturbative Wilson matrix in our program has a geo- 
metric interpretation. The program forms a hollow hypercube about the position of 
the current lattice site for each order. Table 6.1 displays the number of steps from 
the original lattice site (at 0{k)) at which the hypercube is created as a function of k. 
For example, for O(k^) a hollow hypercube of length one is created. The next order 
in K will make a hypercube one step farther out in all directions from the previous 
order, thus, expanding the size of the hypercube by one unit. This is represented in 
the transition between 0{k^) and 0{k^). Our program calculates the scalar loop 

Table 6.1. Dimension of the hollow- hypercubes as a function of k. 





k4 


k5 


k6 




k8 


steps 10 


1 





1 





1 


2 


3 


2 


3* 


2 


3 






4* 


5* 


4 


5 










6 


7 



value to 0{k,^) . According to Table 6.1, the corresponding hollow hypercubes that 
contribute are (1,3,5). For the perturbative VEV calculation the only hypercubes 
that contribute are those that form closed loop, gauge-invariant objects. We see that 
that hypercube of 0{k.^) that has been expanded four steps will not contribute to 
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the sixth order VEV because it is not possible for the gauge-Hnks on this surface to 
form a closed object with links that are three steps away. Therefore, for 0{k.^) no 
hypercube is constructed for the VEV calculation. Similarly, for 0{k^) hypercubes of 
length 3 and 5 are omitted for higher orders in kappa. (All of these values are marked 
with an asterisk.) However, in contrast to the VEV, in the noise calculation these 
contributions are retained because all objects that mimic the noise are included. 

Previous disconnected nucleon calculations have only used a subtraction method 
to reduce the variance of the vector and scalar operators to 0{k.^) and 0{k'^). In this 
thesis, the same calculation is done to 0{nP) and 0{nP), and can easily be extended 
to higher orders in k. 

To determine the perturbative quark propagators we solve the system of equa- 
tions 

MpertX = b, (6.17) 

where b e Z{2). Using a noise vector to solve this system of equations gives quark 
propagators of the off-diagonal (as well as diagonal) elements of MpJ^^. The pertur- 
bative quark propagators from the off-diagonal elements of M^J.^ can have any open 
path up to a given 0{k). These propagators are not gauge invariant and do not con- 
tribute to the operator signal. For example, non-gauge invariant propagator paths of 
0(fi:^) and 0{k,^) are shown in Figure 6.2. 

In practice, we find all the perturbative contributions at each lattice site and 
use that data to reduce the operator variance. 

6.2.2 Vacuum Expectation Values 

The variance-reduction technique reduces the operator signal calculated with 
equation (6.16). To correct for the loss of signal, the perturbative trace must be 
added back into the calculation. 

A picture is instructive to determine which orders in the perturbative expansion 
contribute to the local and vector operators. The scalar and pseudoscalar operators 
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Figure 6.2. Perturbative noise contribution to O(k^) and 0{k'^) used to reduce the variance 
of the disconnected loop operators. 

are local to each individual lattice site. We wish to include all contributions that 
start and end on the same lattice site for these operators. As seen in figure 6.3, the 
orders in kappa which contribute to the local operators are and k,^. These are local 
gauge invariant contributions to the signal. 





4-th Order in Kappa 



6th Order in Kappa 



Figure 6.3. Perturbative scalar operator contributions. 



In general, all even powers of k, contribute to the local operators. 
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For the non-local operators, initial and final lattice sites are connected by a 
gauge-link U{x). In this chapter it is understood that there is an implicit k multi- 
plication in U{x). Thus, in our calculations, contributions to the vector operators 
are of order and k^. In Figure (6.4), the third order diagram is referred to as a 
staple. The fifth order diagram is referred to as a chair diagram. For the vectors, all 
odd-order terms in k contribute. 



3rd Order In Kappm 



51t> Order in Kappa 



Figure 6.4. Perturbative vector operator contribution. 

A similar chair diagram that contributes perturbatively is shown in figure 6.5. 
In this figure, a chair diagram has been constructed around a gauge link U{x). At 
the final position of the object, a "tail" is attached between this lattice site and a site 
that is adjacent. The "tail" is a O(k^) contribution. This contribution is explicitly 
removed by our perturbative subtraction. We refer to this as reducing the variance 
"for free" since there is no extra VEV calculation involved. 

The perturbative vacuum expectation value does not use a noise vector to solve 
the system of equations. Instead we solve the system 

Mx = Ci, (6.18) 
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5th Order In Kappa 

Figure 6.5. Free Order in kappa for vector objects. 

where the vector is an element of the set of all unit vectors that span the Euclidean, 
color, and Dirac spaces. For each Cj, quark propagators are calculated that begin at 
this lattice site and "spread out" to sites that are of 0(k") away. In our calculation 
we determine all the closed loop, gauge-invariant objects up to O(k^) that contribute 
to the operator signal automatically. This is more expensive than an explicit con- 
struction using gauge fields. However, the advantage of this method is that it may 
easily be extended to higher powers in k. 

In our program, there are distinct differences between the noise subtraction part 
and the calculation of the perturbative VEV. The perturbative VEV constructs gauge 
invariant objects that contribute to the signal. To create these gauge invariant objects 
the propagator expands from the current lattice site and only constructs invariant, 
closed objects. The noise subtraction part, on the other hand, constructs all the quark 
propagators to a given order in kappa that contribute to the noise at each site. Hence, 
the difference between the perturbative VEV and the noise-subtraction method is that 
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the VEV propagators are moving outward from a given lattice site to create objects 
that contribute to the signal, while the noise method uses all the contributions from 
every quark propagator ending at the same site. A diagram showing the distinction 
between these two processes is in Figure 6.6. 



Perturbative propagators expand from a given lattice site 
to construct the trace in the vacuum expectation value 
calculation . 



X 



X — * — X 



X 



The perturbative subtraction calculation to reduce the variance 
considers all perturbative propagators that end at a given 
lattice site. 



Figure 6.6. Perturbative VEV and Subtraction Diagrams. 

To create the perturbative VEV contributions mentioned above, a trace over 
Dirac and color indices at each lattice site is used for the perturbative VEV operator 
calculation. The trace is a result of the operator construction in (6.16). In the 
noise-subtraction method, tracing would be incorrect because we wish to determine 
the contribution of the off-diagonal elements to the variance at each lattice site. 
Therefore, in the noise calculation we consider all the quark propagators ending at 
each lattice site as stated above. This information is then used to reduce the variance 
of the exact loop operator. 
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Other perturbative methods have been proposed in references [82,83]. 

6.2.3 Local Operators 

Since the disconnected calculation is delicate, it is always desirable to reduce 
the noise of the operators O. Each calculated operator has a real and imaginary part. 
However, for each operator in question, only the real or imaginary part is needed for 
the calculation. For the Wilson case, operator identities that determine whether the 
real or imaginary part contributes for each operator have been shown using the quark 
propagator identity S — ^58^^^ in reference [6]. The identities are (at each lattice 
site, x): 

Scalar : Re[ip{x)ip{x)] 
Vector : Im['ip{x)jnip{x)] 
Axial : Re[ilj{x)'j5'ynip{x)] (6.19) 
Pseudoscalar : Re[^{x)^^%l){x)\ 
Point — SplitV ector : Klm[^l){x + a^)(l + 'yn)Uj^{x)tl!{x) — tjj{x){l — 7^)?7^(a;)'?/'(a; + a^)] 
Tensor : Im['ip{x)aij,i,ip{x)]. 

In our calculation, the vector and scalar identities play an important role. Although 
these identities are only approximations, when a noise method is employed they al- 
low the omitted part of each operator to be identified as noise and left out of the 
disconnected operator calculation. Ultimately, these identities reduce the variance in 
the form factor calculation. 

6.3 Twisted Mass Disconnected Fermion Loops 
Twisted mass fermions are becoming more popular in hadronic physics because 
they permit calculation of lower quark mass and 0{a) improvement is automatic in 
many quantities [33]. Ultimately, our goal is to use twisted fermions for the nucleon 
and disconnected operators to make realistic calculations of nucleon strange form 
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factors. This section presents a method in which the perturbative subtraction method 
can be implemented to determine disconnected loop values using twisted fermions. 

To make use of noise methods for twisted mass fermions we need to find an 
equivalent matrix to that in equation (6.13). The twisted mass ferimonic matrix is 

Mtm = I - 2ifx-f5T3 - kP, (6.20) 

where the matrix P is the same matrix used in the Wilson formulation. Let tan(5) = 
2kij,. The twisted fermion matrix can then be written as 

Mtm = / - i75 tan(5) - kP. (6.21) 

The inverse of this matrix can be shown to be 

M^Ij = e'^^' cos{5){- „ ). (6.22) 

In this form we can expand equation (6.22) in a geometric series such that perturbative 
subtraction may be employed. The expansion of the twisted fermion matrix is 

Mtm{P) = cos(5)e*^^5(/ + cos(5)Ke*^^^P + cos2(5)K2e*^^5Pe'^^s (6.23) 
+ cos^ (5) K^e'^'^' Pe**^« Pe^*^« P + ...). 

The extension from the Wilson perturbative subtraction method to the twisted 
fermion method amounts to a chiral rotation of the original quark matrix. This 
expansion of the quark matrix is referred to as a left twisted mass expansion because 
the chiral twist projects onto the Wilson matrix from the left hand side. An equivalent 
expression can be formed by expanding the quark matrix such that the the rotation is 
applied from the right. These are numerically identical, but the right multiplication 
was found to be more expensive in terms of computational time. The extra time is a 
result of the rotation being done after each hypercube has been formed instead of a 
simple chiral rotation on the initial noise vector. 
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The quark charge conjugation property in the full twisted mass formalism is 
Du = 75^)^75 [39]. Using this property we can show the local operators to have 
the same form as in equation (6.19) when averaged over "tmU" and "tmD" quarks. 
The calculation is different in that one has two flavors of quark propagators. In 
this formalism the charge conjugation property changes flavor as well as charge. For 
example, the scalar operator is only completely real when both flavors are considered 
in the trace, 

O^^ = ^ < ipuip + ipdipd > (6.24) 
= l(M-^ + M-i). (6.25) 

It is important to realize that the {u, d) subscripts on the quark propagators are 
not the physical up and down quark. Instead they are the unphysical twisted mass 
labels [84]. With this in mind, we now have two flavor doublets on the lattice. The 
doublets can be written as 

i^i^Q,^k^C-), (6.26) 

where the subscripts {I, h) are for light and heavy respectively. Each doublet is mass 
degenerate, thus the c quark is not the physical charm quark. Since this "charmed 
quark" is not used in any of our quenched lattice calculations it is not necessary to 
include an explicit nondegeneracy in the doublet. This procedure is employed by the 
authors of [85] . 

6.4 Subtraction Results 
The first results of higher order subtraction in the twisted mass basis are pre- 
sented in this section. In the Wilson case, for much heavier quark masses = .148), 
it has been shown that the effects of higher order subtraction can be dramatic for 
point-split vector, vector, and scalar operators [6]. The estimated computer time 
used to do a perturbative subtraction calculation of the disconnected operators is 
determined by the ratios of unsubtracted variance to the subtracted variance. 
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Table 6.2. Sixth Order k Variance Ratios of Scalar Operators. 



Subtraction 


Scalar 


PseudoScalar 


F ourthOrder{K^ ) 


1.9 


2.1 


SixthOrder{K^) 


2.0 


2.1 



The ratios of the variance were calculated for k'^ and with fifty twisted 
mass configurations to investigate the computational gains from the higher order 
subtraction. The hopping and twisted mass parameter for this calculation are fj, = 
0.30 and k = 0.15679, respectively. Each configuration uses the optimum number of 
noises in the determination of the perturbative quark matrix M^l^{P). 

The optimum number of noises to minimize the variance can be determined with 
the variances of the gauge configurations and noises, Vgauge and Vnmsei respectively [6]. 
Given A^-configurations and M-noises per configuration, the error bar on a given 
operator is 

V NM N ^ ^ 

Clearly, equation 6.27 is minimized for M = 1. This result can be modified to incor- 
porate computational overhead. If it is assumed that there is an overhead associated 
with generating configuration and we assume a fixed amount of computer time for 
each configuration, then 

T = NM + GnN, (6.28) 

where is the time overhead for configuration generation. The minimization of 
equation (6.28) gives 

M = ^^y^^ (6.29) 

'-Vgauge 

where Snoise and Sgauge are determined by their respective variances. In our calculation 
the ratio Snoise/ Sgauge ~ 1 for the vectors, which results in an optimum number of 
noises of approximately 5. The number of noises was not optimized for the scalar but 
the vector operators. 
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Table 6.3. Sixth Order k Variance Ratios of Vector Operators. 



Subtraction 


ChargeDensity 


Jlcurrent 


J2current 


J3current 


FourthOrderRatio{K^) 


5.0 


4.3 


3.8 


4.0 


SixthOrderRatio{K^ ) 


5.9 


5.0 


4.6 


4.5 



In table 6.2, it is observed that the scalar operators do not respond as well to 
the perturbative subtraction as do the currents. This behavior is consistent with what 
is found in the Wilson case. The scalar operators gain a factor of 1.9 in computer 
time using fourth order subtraction. In comparison, when sixth order subtraction 
is used for the scalar operator a gain of approximately 2.0 is reported. This is an 
approximate 10 percent increase in algorithm speed. 

According to table 6.3, using higher order subtraction one saves a factor of 
approximately five in computer time for the spatial currents. An equivalent statement 
is that the number of noises needed to produce a comparable result to unsubtracted 
noise method is reduced by a factor of five. The charge density responds better to the 
subtraction method than the spatial currents. The charge density operator saves a 
factor of approximately six in computer time. These results indicate an approximate 
20 percent increase in algorithm speed from 0{k'^) to 0{k^). 

The scalar operator signal in figure 6.7 increases at time steps 1 and 32. This 
edge effect is a result of the non-periodic boundary condition in the time direction. 
Fortunately, these values are not used in the correlation function calculation and can 
be ignored. 

Similar subtraction diagrams for the psuedoscalar(Figure 6.8), charge density 
(Figure 6.9), and a spatial current (Figure 6.10) are below. These figures support the 
conclusion that the twisted vector operators respond better to subtraction methods 
than the scalar operators. 

In the tmQCD formalism it has been shown that the scalar-pseudoscalar and 
axial vector- vector operators mix [13]. In this thesis the nucleons are calculated at 
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Figure 6.7. Scalar Subtraction to 6th order in k. 

maximal twist in which the physical constraint that there is no mixing between the 
charged psuedoscalar and vector is imposed to eliminate the axial- vector mixing [23]. 
However we have found that this does not eliminate the scalar-psuedoscalar operator 
mixing. Scalar mixing was observed by the authors of reference [23]. This mixing can 
be seen in equation 6.23. The first term in the expansion is approximately 1 for a small 
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Figure 6.8. PseudoScalar Subtraction to 6th order in k. 
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Figure 6.9. Charge Density Subtraction to 6th order in n. 




Figure 6.10. Ji-current Subtraction. 

rotation angle S. This guarantees that the scalar operator has a vacuum expectation 
value. So, independent of the maximal twist angle the scalar-pseudoscalar mixing will 
occur with this approach. Hopefully, other methods can be determined to remove the 
scalar-pseudoscalar mixing and promote a twisted disconnected noise method. 

Using subtraction methods give a useful tool for exploring quark loops in the 
disconnected sector. In future disconnected nucleon calculations hopefully twisted 
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perturbative subtraction can be used where scalar-pseudoscalar mixing has been re- 
moved to reach hghter quark masses. 



CHAPTER SEVEN 
Numerical Simulations and Results 

The gauge field configurations used in this study were generated from the unim- 
proved Wilson gauge action at /3 = 6.0 on a 20^ x 32 lattice corresponding to a lattice 
spacing of 

a = 0.1011(7)/m. (7.1) 

as obtained from reference [86] from a physical string tension of \fK =427MeV. This 
lattice spacing was used in the strangeness calculation in [1]. Our full ensemble of 
200 configurations was produced from a thermalized Markov chain. Each ensemble 
configuration is generated with 2000 heatbath updates between saved configurations. 

The twisted mass lattice action at maximal twist was used to obtain four valence 
quark masses per configuration. Each mass has an associated hopping parameter k 
and twisted mass parameter \x. These values are in Table (7.1). 



Table 7.1. Maximally twisted mass pairs, (k,/x). 



Mass Number 


Hopping parameter, k 


Twisted mass parameter, 


1 


0.15679 


0.030 


2 


0.15708 


0.015 


3 


0.15721 


0.010 


4 


0.15728 


0.005 



The valence quarks in our simulation are subject to Dirichlet time boundaries. 
The source is located at (1,1,1,4) which is four timesteps away from the boundary. 

Strangeness matrix elements are calculated using standard methods in which 
the three-point function is created by correlating a strange-quark loop with the nu- 
cleon propagators. The strange-quark loops are calculated with the perturbative 
subtraction techniques from chapter 6 with real Z-i noise. The scalar loops in our 
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calculation are determined to O(k^). The previous nucleon strangeness calculation 
employing this stochastic technique (Ref. [1]) was to lowest order subtraction, 0{k^), 
from reference [87]. 

In our disconnected calculation we use loop values of Kioop = 0.152 and Kioop = 
0.154. These kappa values correspond to vector meson masses of 912(8) MeV and 
1066(4) MeV respectively [1, 86] which surrounds 1019 MeV which ensures that our 
data will interpolate to a strange quark loop. In our lattice simulation the matrix 
elements are extracted from the ratio in equation (4.36). A fixed loop background 
starting at the source and ending at time step 20 was used in these calculations. We 
consider the lowest five momentums given by, 

aV = n(7r/10)^n = 0,l,2,3,4. (7.2) 

The lowest three momentum are focused on in this chapter because the momentum 
associated with n = 3, 4 are still noisy and unpredictable due to a lack of configura- 
tions. 

Figures 7.1 - 7.8 plot the ratio of three to two-point functions for the lowest 
three momentums of all four masses in table (7.1). 

7. 1 Jackknife Error Bars and Linear Fit 
The uncertainties for all the figures in this chapter are calculated with a jack- 
knife error bar technique [88,89]. A jackknife error bar is calculated for the ratio 
in equation (4.36) at every time step. The jackknife technique is referred to as a 
resampling method because it uses small changes from the original data set to deter- 
mine the uncertainty in the data. The change in the data is a result of omitting each 
configuration from the ensemble average one-by-one and reproducing the ratio of the 
three and two-point functions at each time slice. The jackknife method is summarized 
below. If we define the jackknife averages of the ratios from equation (4.36) to be the 
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Figure 7.1. Ratio in equation (4.36) for the first three momenta for mass 1 scalar density 
diagram at Kioop=-^52. 
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Figure 7.2. Ratio in equation (4.36) for the first three momenta for mass 1 scalar density 
diagram at K;oop=- 154. 



configuration average of the ratios while omitting the z -ratio, then, 
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Figure 7.3. Ratio in equation (4.36) for the first three momenta for mass 2 scalar density 
diagram at Kioop=-^52. 
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Figure 7.4. Ratio in equation (4.36) for the first three momenta for mass 2 scalar density 
diagram at K;oop=- 154. 

Also, we define the jackknife estimate of the ratio as the configuration average over 

all jackknife averages defined in equation (7.3), 

1 ^ 

i=l 
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Figure 7.5. Ratio in equation (4.36) for the first three momenta for mass 3 scalar density 
diagram at Kioop=-^52. 
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Figure 7.6. Ratio in equation (4.36) for the first three momenta for mass 3 scalar density 
diagram at K;oop=- 154. 



The uncertainty in the ratio is then, 

a^J = 7^(< {R'^y>-{< R' >n (7.5) 
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Figure 7.7. Ratio in equation (4.36) for the first three momenta for mass 4 scalar density 
diagram at Kioop=-^52. 
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Figure 7.8. Ratio in equation (4.36) for the first three momenta for mass 4 scalar density 
diagram at K;oop=- 154. 

The values of a^j are the uncertainties calculated for each time-slice and momentum 
for all masses in the figures above. 
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Linear fits of the data in figures (7.1) - (7.8) are needed to extract the matrix 
elements. A least-squares fit to an arbitrary function is used to find the best fit over 
a specific range of time-steps [90]. Let x — xi^^^) be a fit parameter that is a 
function of the uncertainties generated with the jackknife technique and the ratio, R, 
from equation (4.36). The method then searches the parameter space of x to fiiid 
its minimum value which corresponds to the best uncorrelated linear fit of the data. 
This method was adapted to consider correlated fits of the jackknife ratio data [91]. 
The fit parameter is multiplied by the covariance matrix, Cij, defined in equation 



where N is the number of configurations, represent different time slices, Rfin) 

has jackknife ratio data, and < > is the configuration average that removes the 
bias. The covariance matrix considers correlations between ratio data at different 
time slices. These correlated fits are used in this thesis because it predicts the best 
linear fit over a specific time interval and the corresponding uncertainty in that fit. 
The fits and error bars presented in Table (7.2) and Table (7.3) are calculated with 
this correlated least-squares method. 



The data in Tables (7.2) and (7.3) are for K/o„p = .152 and Kioop = .154, respec- 
tively. In each table, the lowest three momentums are reported for all four twisted 
mass {k, /i) pairs. The range of the time interval for the linear least-squares fit and 
the associated error bar are given. The best hnear fits for each mass and momentum 
are reported. Since this is a low statistics study giving preliminary results, the range 
of each fit is different for each mass and momentum. 

The value of the scalar for the hghtest quark mass at second momentum for 
Kioop — 0.154 is considered over the 15 — 19 time slices in Table (7.3). Figure (7.8) 



(7.6). 




(7.6) 



n=l 



7.2 Discussion 
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Table 7.2. Fits for the matrix elements from equation (4.36) for all 4 masses 
and Kioop = 0.152. The momentum squared is a^g^ = n(7r/10)^. 





time steps 


n 


Scalar 


(0.15679, 0.030) 










13-17 





2.1 ± 1.5 




15-17 


1 


0.77 ± 0.64 




15-19 


2 


0.37 ± 0.30 


(0.15708, 0.015) 










13-16 





2.3 ± 1.9 




15-18 


1 


1.2 ± 0.89 




12-16 


2 


0.88 ± 0.81 


(0.15721, 0.010) 










15-19 





2.5 ± 1.8 




14-17 


1 


1.6 ± 1.2 




9-13 


2 


0.90 ± 0.76 


(0.15728, 0.005) 










13-18 





3.2 ± 3.0 




14-18 


1 


1.3 ± 1.2 




9-12 


2 


0.97 ± 0.79 



Table 7.3. Fits for the matrix elements from equation (4.36)for all 4 masses 
and Kioop — 0.154. The momentum squared is — n(7r/10)^. 



{Ky,iJ,) time steps n Scalar 

(0.15679, 0.030) 



(0.15708, 0.015) 



(0.15721, 0.010) 



(0.15728, 0.005) 



13-16 





14-18 


1 


15-19 


2 


11-16 





14-18 


1 


14-18 


2 


11-16 





14-18 


1 


13-19 


2 


9-13 





15-19 


1 


10-13 


2 



2.1 ± 1.6 
1.0 ± 0.75 



0.49 


± 


0, 


.38 


1.9 


± 


1, 


.5 


1.4 


± 


1, 


.0 


0.92 


± 


0, 


.74 


2.1 


± 


1, 


.7 


1.8 


± 


1, 


,4 


1.6 


± 


1, 


.3 


2.3 


± 


2, 


.0 


1.3 


± 


1, 


.0 


0.97 


± 


0, 


.89 



shows the plot of this data. The data point at the 19 time slice is included because 
the fitting routine suggests that this point is highly correlated over this range and 
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reasonable to fit with. The scalar value, without this point, in the range 15 — 18 is 
3.4 ± 3.3. This large change in the scalar value is a result on insufficient statistics 
and will be resolved with the addition of more configurations. 

Plots of the scalar ipip as a function of the dimensionless 4-momentum transfer 
squared (a^Q^) are plotted in Figures (7.9) and (7.10). The square of the 4-momentum 
transfer is 

Q' = {q- q)\ (7.7) 

where q = [E, q) and q' = {rriN, 0, 0, 0) is the final and initial momentum respectively, 
and m-AT is the nucleon mass from [23]. Then can be written 

= 2m{E -m). (7.8) 

The nucleon masses used in the 4-momentum transfer plots were calculated in ref- 
erence [23]. One can see that the scalar density falls of smoothly and has similar 
behavior for both disconnected loop values. 




Figure 7.9. The scalar density for the first three momentums and all momentums as a 
function of the dimensionless 4-momentum transfer a?Q'^. These plots are for Kioop=-^52. 
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Figure 7.10. The scalar density for the first three momentums and all momentums as a 
function of the dimensionless 4-momentum transfer a^Q^. These plots are for i^ioop — • 

154. 

The zero momentum values in Figures (7.9) and (7.10) are offset slightly so that 
the data points can be clearly identified. 

Another useful plot is the scalar density as a function of the pion mass. The 
pion mass squared is proportional to the quark mass. The pion mass is used in xPT 
to extrapolate to the physical quark masses. A plot of the scalar {ipip) as a function 
of the pion mass squared for each of the loop values is found in Figure (7.11) and 
Figure (7.12). The pion masses corresponding to our maximal twist masses used in 
these plots are reported in [23]. 

A similar calculation of the strangeness scalar density matrix elements in the 
Wilson formalism is presented in reference [1]. Their calculation was performed on 
a 20^ X 32 lattice with (3 = 6.0. The valence quark hopping parameters were = 
{0.152, 0.153, 0.154} with loop values of nioop = 0.152 and Kioop = 0.154. The results in 
this high statistics calculation are determined with 2000 configurations with statistical 
uncertainties obtained from 3000 bootstrap ensembles. The fit to the scalar density 
in this paper are found in Table (7.4). The fits in this table begin 10 time steps from 
the source. 
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Figure 7.11. The zero momentum scalar density as a function of pion mass squared. These 
plots are for Kioop=-^52. 



Figure 7.12. The zero momentum scalar density as a function of pion mass squared. These 
plots are for K;/oop=-154. 



The scalar results in [1] were found to decrease in amplitude as the momentum is 
increased. The scalar values in Tables (7.2) and (7.3) also decrease as the momentum 
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Table 7.4. Fits for the scalar matrix elements from reference [1]. 
These scalar values are calculated for Kioop— 0-152 and Kioop—O.lbA 





n 


Kioop = 0.152 


Kioop = 0.154 


0.152 













2.6(4) 


3.7(13) 




1 


1.7(2) 


2.1(6) 




2 


1.2(2) 


1.1(6) 


0.153 













2.7(5) 


4.0(14) 




1 


1.8(3) 


2.2(7) 




2 


1.3(2) 


1.3(11) 


0.154 













2.9(5) 


4.2(5) 




1 


1.8(3) 


2.3(8) 




2 


1.3(3) 


1.3(8) 



increases. The scalar data reported in the high statistics study was found to be 
independent of the valence quark masses. Our data appears to increase slightly for 
smaller valence quark mass for both loop values. The nucleon quark masses in our 
study are the lightest valence masses used for the nucleon strangeness calculation to 
date and, therefore, our preliminary results are the first twisted mass calculation of 
the nucleon strangeness scalar density. 

The lightest valence quark mass in [1] is most comparable to the heaviest mass 
in our simulation. The matrix elements amplitudes of the lowest three momenta 
for the twisted and Wilson case are different. For example, in the Wilson case, the 
zero momentum scalar densities for — 0.154 with Kioop — 0.152 and Kioop — 0.154 
were reported to be 2.9 and 4.2 respectively. Our scalar values most comparable to 
Kv = 0.154 arc 2.1 for both loop values. Our heaviest twisted quark mass pair is 
{k — 0.15679, /I — 0.030). The values for the higher momentums compare similarly. 

Our preliminary results suggest that the raw data for the scalar elements are 
being calculated correctly. This calculation is aimed toward forming the renormaliza- 
tion group invariant quantity representing the fractional strange quark contribution 
to the nucleon mass in equation (7.9). 
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rris < N\ss\N > (0) 

Once higher statistics are acquired the physical masses can be obtained in the con- 
tinuum hmit using xPT. The success of this calculation gives hope for future, high 
statistics, calculations using these methods of the electric and magnetic form factors 
to determine electric and magnetic properties of the nucleon in the presence of a 
strange quark loop. 



CHAPTER EIGHT 
Conclusion 

A study of the strangeness contribution to the nucleon was conducted in this 
thesis. Our results show that the methods presented here are viable and will allow 
for a better study of the strangeness content of baryons. More specifically, we have 
shown preliminary results that indicate that the twisted mass formalism is a good 
approach to calculate the scalar form factor. To calculate the scalar many new and 
interesting techniques were developed to "zero-in" on the form factor values using 
lighter valence quark masses so that we can make better contact with experimental 
results. Future work will include a high statistics calculation of the electric and 
magnetic form factors so that one may have a better understanding of the nucleon 
electromagnetic properties. 

We have shown many techniques to improve lattice calculations. Our results 
show that useful variations of the GMRES(m) algorithm can be employed to solve 
systems of linear equations that arise in Lattice QCD calculations efficiently. The 
saved matrix-vector products from these algorithms can reduce computational time 
dramatically over the life of a high statistics lattice calculation. Specifically, we 
have shown that GMRES-DRS(m,k) is a good technique to solve shifted systems of 
equations in the Wilson case by taking advantage of the properties of the Krylov 
subspace. As an extension to GMRES-DRS(m,k), we have developed another new 
technique to use a shifted GMRES(m)-Proj(k) method to solve subsequent right-hand 
simultaneously after the base system has used GMRES-DRS(m,k). 

The disconnected quark loop calculation used to form the disconnected three- 
point function was improved by expanding to higher orders in the perturbative ex- 
pansion. This is an important result because going to higher order in kappa fur- 
ther reduces the variance of the loop operators and saves valuable computer time 
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in the calculation. A twisted mass noise method was also presented in this thesis. 
This method responds well to the subtraction techniques in that the variance of the 
twisted loop operators is significantly reduced in our simulations. As noted in chapter 
6, these loops suffer from scalar-pseudoscalar mixing that causes both the scalar and 
the pseudoscalar to acquire a VEV. Future work to remove the mixing in the twisted 
perturbative subtraction method is necessary so that one may go to lower quark mass 
for the loops and produce more accurate strangeness calculations. 
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